Métodos Numéricos para Resolu¸cão de Equa¸cões ... · Resumo Neste trabalho desenlvolvemos um método de resolu¸cão de problemas de valor inicial ... Um artigo recente,

UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE

CENTRO DE CIENCIAS EXATAS E DA TERRA

PROGRAMA DE POS-GRADUACAO EM MATEMATICA

APLICADA E ESTATISTICA

Metodos Numericos para Resolucao de Equacoes

Diferenciais Ordinarias Lineares Baseados em

Interpolacao por Spline

Thiago Jefferson de Araujo

Orientador: Prof. Dr. Nir Cohen

Natal, agosto de 2012

UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE

CENTRO DE CIENCIAS EXATAS E DA TERRA

PROGRAMA DE POS-GRADUACAO EM MATEMATICA

APLICADA E ESTATISTICA




Thiago Jefferson de Araujo

Dissertacao de mestrado apresentada ao Programa

de Pos-Graduacao em Matematica Aplicada e Es-

tatıstica da Universidade Federal do Rio Grande

do Norte (PPGMAE - UFRN) como requisito par-

cial para obtencao do tıtulo de Mestre em Ma-

tematica Aplicada e Estatıstica.





Este exemplar corresponde a redacao final, defen-

dida por Thiago Jefferson de Araujo e aprovada

pela comicao julgadora.


Banca Examinadora:

• Prof. Dr. Nir Cohen

• Prof. Dr. Allan de Medeiros Martins

• Prof. Dr. Daniel Cordeiro de Morais Filho

Resumo

Neste trabalho desenlvolvemos um metodo de resolucao de problemas de valor inicial

com equacoes diferenciais ordinarias baseado em splines, com enfase em equacoes lineares.

O metodo serve como alternativa para os metodos tradicionais como Runge-Kutta e no

caso linear com coeficientes constantes, evita o calculo de raızes de polinomios. O metodo

foi aplicado para um problema central da teoria de controle, o problema de resposta a

degrau para uma EDO linear, incluindo o caso de coeficientes nao-constantes, onde a al-

ternativa pelo calculo de raızes nao existe. Implementamos um algoritmo eficiente que usa

apenas operacoes tipo matriz-vetor. O intervalo de trabalho (ate o tempo de acomodacao)

para as equacoes estaveis com coeficientes constantes e determinado pelo calculo da raiz

menos estavel do sistema, a partir de uma adaptacao do metodo da potencia. Atraves de

simulacoes, comparamos algumas variantes do metodo. Em problemas lineares gerais com

malha suficientemente fina, o novo metodo mostra melhores resultados em comparacao

com o metodo de Euler. No caso de coeficientes constantes, onde existe a alternativa ba-

seada em calculo das raızes, temos indicacoes que o novo metodo pode ficar competitivo

para equacoes de grau bastante alto.

PALAVRAS CHAVE: Equacao Diferencial Ordinaria Linear; Estabilidade Assintotica

e Transiente; Metodo da Potencia; Metodos de Interpolacao por Spline.

i

Abstract

In this work we have elaborated a spline-based method of solution of inicial value

problems involving ordinary differential equations, with emphasis on linear equations.

The method can be seen as an alternative for the traditional solvers such as Runge-Kutta,

and avoids root calculations in the linear time invariant case.

The method is then applied on a central problem of control theory, namely, the step

response problem for linear EDOs with possibly varying coefficients, where root calcula-

tions do not apply. We have implemented an efficient algorithm which uses exclusively

matrix-vector operations. The working interval (till the settling time) was determined

through a calculation of the least stable mode using a modified power method.

Several variants of the method have been compared by simulation. For general linear

problems with fine grid, the proposed method compares favorably with the Euler method.

In the time invariant case, where the alternative is root calculation, we have indications

that the proposed method is competitive for equations of sifficiently high order.

KEY WORDS: Linear Ordinary Differential Equation; Asymptotic Stability and

Transient; Power Method; Spline Methods.

ii

Dedicatoria

”Esse trabalho e dedicado a minha mae pelo amor

incondicional e por ter sempre me motivado a es-

tudar.”

iii

Agradecimentos

Agradeco ao professor Nir Cohen (UFRN) pela orientacao, pelos concelhos pessoais e

por compartilhar um pouco da genialidade e humor.

Agradeco a Anna Rafaella por estar ao meu lado durante todo o perıodo que fui aluno

do PPGMAE, pelo carinho e por todo apoio que me deu me ajudando sempre que precisei.

Ao professor Marcelo Gomes (UFRN) por tudo que me ensinou durante a minha

graduacao, pelo exemplo de professor e pessoa que eu tanto tento me espelhar e por

entender as dificuldades que passei na epoca, muito obrigado.

A professora Viviane Simioli (UFRN) por ter me ajudado em diversos momentos alem

do perıodo em que fui seu aluno e por ser uma excelente professora.

Ao professor Roberto Hugo (UFRN) que durante a disciplina de Otimizacao me levou

a estudar assuntos que foram fundamentais para que eu pudesse escrever esse trabalho.

Ao professor Allan de Medeiros (UFRN) por participar da minha qualificacao, da

minha defesa e por ter ajudado diversas vezes na construcao desse trabalho.

A professora Debora Borges (UFRN) por ter participado da minha banca qualificacao.

Ao professor Daniel Cordeiro (UFCG) por participar da banca da minha Defesa e

pelas sugestoes.

A CAPES por ter me concedido bolsa durante boa parte do mestrado.

iv

Aos meus amigos que foram membros do PET com quem eu pude aprender e me

divertir muito, em particular Jobson Hugo e Ailton Rodrigues.

Aos servidores do CCET.

A todos os meus amigos e familiares.

v

Lista de Tabelas

2.1 Iteracoes do metodo da Potencia aplicado a matriz do exemplo 2.1 . . . . . 31

2.2 Iteracoes do metodo da Potencia Adaptado aplicado a matriz do exemplo 2.2 38

2.3 EDO’s lineares com coeficientes constantes para estudar o tempo de aco-

modacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.1 Resolucao do exemplo 3.5 pelo metodo da Suavidade com s = 100 e xf = 20 69

3.2 Comparacao entre as solucoes do exemplo 3.7 pelo metodo de Euler e pelo

metodo da Suavidade com o calculo dos autovalores . . . . . . . . . . . . . 77


metodo da Suavidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78


metodo da Suavidade e pelo calculo dos autovalores . . . . . . . . . . . . . 79


metodo da Suavidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

vi

Lista de Figuras

1.1 Grafico da solucao do exemplo 1.5 . . . . . . . . . . . . . . . . . . . . . . . 14

1.2 EDO (exemplo 1.6) onde a parte real dos autovalores e nula . . . . . . . . 17

1.3 EDOs estaveis (a parte real de todos os autovalores e negativa) . . . . . . . 18

1.4 EDOs instaveis (a parte real de um dos autovalores e positiva) . . . . . . . 19

1.5 Solucao das EDOs das equacoes 1.14, 1.14, 1.16 e 1.17 . . . . . . . . . . . . 22

2.1 Analise do tempo de acomodacao das EDO’s da tabela 2.4 . . . . . . . . . 43

3.1 Grafico da solucao do exemplo 3.1 . . . . . . . . . . . . . . . . . . . . . . . 48

3.2 Spline cubico para aproximacao do calculo de uma EDO . . . . . . . . . . 53

3.3 Solucao Exata (tracejado) e spline calculado pelo metodo da Sub-particao . 56

3.4 Solucao Exata (tracejado) e spline calculado pelo metodo da Suavidade . . 61

3.5 Solucao Exata (tracejado) e spline calculado pelo metodo da Suavidade

Implementado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

3.6 Exemplo de como aplicar o metodo da Suavidade Adaptado a uma EDO

de grau 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

vii

3.7 Solucao Exata (tracejado) e spline calculado pelometodo da Suavidade

Adaptado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

3.8 Solucao da EDO do exemplo 3.7 pelo metodo de Euler com s = 175 e da

Suavidade com s = 70 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

3.9 Solucao da EDO do exemplo 3.8 pelo metodo de Euler e da Suavidade . . . 80

viii

Sumario

Introducao 2

1 Equacoes Diferenciais Ordinarias Lineares 5

1.1 Transformada de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2 Fracoes Parciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3 Solucao e Analise da Estabilidade de EDO’s Lineares Com Coeficientes

Constantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4 Pontos Crıticos da Solucao de EDO’s Lineares com Coeficientes Constantes

de Grau 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

1.5 Sistemas de Equacoes Lineares de Primeira Ordem . . . . . . . . . . . . . 23

2 Metodo da Potencia 25

2.1 Metodo da Potencia para Matrizes com um Unico Autovalor Dominante . . 27

2.2 Metodo da Potencia para Matrizes com um Par de Autovalores Dominantes

Complexos Conjugados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.3 Metodo da Potencia para o Calculo do Autovalor comMaior Valor Absoluto

de uma Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

ix

2.4 Tempo de Acomodacao de EDO’s Lineares Com Coeficientes Constantes . 40

3 Metodos Splines para o Calculo Aproximado da Solucao de EDOs 45

3.1 Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.2 Metodo da Sub-particao . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3 Metodo da Suavidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.4 metodo da Suavidade Adaptado . . . . . . . . . . . . . . . . . . . . . . . . 70

3.5 Resultados Numericos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Conclusao 82

A Propriedades Basicas da Transformada de Laplace 84

B Estabilidade Assintotica de EDOL’s com Coeficientes Constantes 86

C Programacao no MATLAB 88

Referencias Bibliograficas 94

Introducao

Um problema de valor inicial com Equacao Diferencial Ordinaria de Linear com Co-

eficientes Constantes pode ser resolvido pelo Metodo de Euler (ou metodos gerais de

Runge-Kutta), ou a partir do conhecimento do polinomio caracterıstico da equacao (ver,

por exemplo, Boyce e Di Prima [2] e o capıtulo 1 desse trabalho). Porem esse ultimo

metodo e algo limitado a EDO’s Lineares com Coeficientes Constantes e mesmo nesse

caso o custo computacional para o calculo dessas raızes pode ser alto, a depender do grau

da equacao.

Neste trabalho sera desenvolvido um metodo alternativo, baseado na teoria de splines,

que evita o calculo das raızes e e numericamente competitivo com o metodo de Euler.

A motivacao deste trabalho e a analise de resposta ao degrau, e em particular do regime

transitorio, dentro da teoria de sistemas dinamicos na engenharia, usado na analise de

funcoes de transferencia para sistemas lineares com coeficientes constantes. Isto e baseado

no fato que a resposta a degrau define completamente a funcao de transferencia (Lathi

[7]), ou seja, permite recuperar completamente o sistema linear.

O estudo de sistemas estaveis trata, de um lado, o regime transitorio, que engloba

a solucao em tempo limitado, e o regime permanente, onde o sistema aproxima-se do

equilıbrio, ou seja, do seu valor assintotico. Ressaltamos que a analise do regime tran-

sitorio e uma parte importante de algumas aplicacoes, por exemplo, em sistemas de alta

tensao, agendamento de controladores (gain scheduling) na aviacao, ou o estudo de choque

e amortecimento para sistemas de varios tipos.

2

3

Em sistemas lineares, enquanto a parte real dos autovalores define o regime permanente

e a estabilidade assintotica, o regime transitorio e pouco estudado em termos espectrais

e normalmente e calculado numericamente. Um artigo recente, Cohen e Lewkowicz [4],

liga a qualidade do regime transitorio a riqueza do conjunto de funcoes de Lyapunov do

sistema, embora ainda falte testes numericos para essa afirmacao.

Neste trabalho adotamos uma abordagem diferente para a resolucao de problemas

de valor inicial, baseada em splines polinomiais. Esta abordagem permite a descricao

numerica eficiente da solucao no regime transitorio, para sistemas estaveis, e com isso,

uma analise qualitativa da solucao. Apesar de ser um metodo geral, assim como o metodo

de Euler, capaz de tratar equacoes gerais, o metodo e particularmente bem comportado

no caso de equacoes lineares invariantes no tempo (LIT).

Splines sao usados principalmente em problemas de interpolacao, onde uma curva

suave deve ser encontrada, passando em alguns pontos dados. O spline e determinado

atraves de dois tipos de restricoes: as restricoes de interpolacao nos pontos dados e res-

tricoes de suavidade entre duas secoes vizinhas do spline. Aqui, no lugar de condicoes

de interpolacao, as condicoes forcam a curva a satisfazer a EDO em alguns pontos da

curva. Na medida em que pudemos verificar, nao encontramos esta ideia natural do uso

de splines na literatura.

O metodo de splines requer a solucao de sistema linear em cada passo, para determinar

os coeficientes de um k-esimo polinomio. Aqui, apresentamos e comparamos duas vari-

antes que chamamos o Metodo da Sub-particao e o Metodo da Suavidade. Comparando

os dois, o metodo de Sub-particao apresenta menos suavidade, mas melhor aderencia a

equacao diferencial.

Da nossa experiencia (veja Capıtulo 3), o metodo de suavidade e numericamente me-

lhor que o da Sub-particao e, ao mesmo tempo, permite implementacao mais simples.

Especificamente, a implementacao numerica dos dois algoritmos requer apenas multi-

plicacao matriz-vetor (Golub e Van Loan [5]) no caso de EDO LIT; e a matriz e a mesma

para todo k.

4

Nosso metodo evita completamente a necessidade de calcular raızes de um polinomio.

Porem, uma opcao viavel e a analise do modo principal de estabilidade, ou seja, a raiz

com a maior parte real. Com isto, pode-se ser calculado o tempo de acomodacao do

sistema, onde o sistema passa para o regime permanente, permitindo a determinacao do

intervalo de trabalho para o metodo de splines. Aqui existem duas possibilidades: caso o

modo principal for real, ele e calculado usando o metodo da potencia aplicado na matriz

associada com a EDO (ver, por exemplo, Poole [12]). O caso de duas raızes complexas

conjugadas requer uma adaptacao do metodo da potencia, que apresentamos em Capıtulo

2.

No ultimo capıtulo incluımos alguns testes numericos que estudam o tempo computa-

cional e erro relativo: uma comparacao entre os Metodos da Suavidade e de Euler para

EDOs Lineares; e no caso LIT, uma comparacao com o metodo explıcito baseado no

calculo das raızes.

Capıtulo 1

Equacoes Diferenciais Ordinarias

Lineares

1.1 Transformada de Laplace

Nesse trabalho a Transformada de Laplace sera usada unicamente com o intuito de

resolver Equacoes Diferenciais Ordinarias Lineares, visto que ela permite resolver essas

equacoes a partir de operacoes com equacoes polinomiais que e bem mais simples de se

realizar, porem ela pode ser usada em diversas situacoes. Pode encontrar mais detalhes

sobre o assunto Ogata [10] e em Madureira [8], trataremos aqui apenas os resultados

essenciais para o desenvolvimento do trabalho.

Definicao 1.1. Se y e uma funcao complexa, definida para todo t ≥ 0 e s e um numero

complexo onde s = p+ qi, com p, q ∈ R donde para cada p > 0 a integral

Y (s) = Ly(s) = L[y(t)] =

∫ ∞

0

e−sty(t)dt (1.1)

e convergente, entao a funcao Y e chamada de transformada de Laplace da funcao y.

A Transformada de Laplace de uma funcao satisfaz as seguintes propriedades:

5

Capıtulo 1. Equacoes Diferenciais Ordinarias Lineares 6

p1 : L[a1y1 + a2y2] = a1L[y1] + a2L[y2], (linearidade)

p2 : L[y(k)] = skL[y]− sk−1y(0)− sk−2y(1)(0)− ...− sy(k−2)(0)− y(k−1)(0)

p3 : L[eat] =

1

s− a

onde a, a1, a2 ∈ C e y, y1 e y2 sao funcoes definidas para todo t ≥ 0 e y(i) denota a

i-esima derivada de y.

Para calcular o valor de y(t) = eλt quando o valor de λ e complexo, sera aplicada a

Relacao de Euler.

Relacao de Euler: Para todo numero real r vale a seguinte relacao

eri = cos(r) + i sen(r). (1.2)

Se λ = c+ di, onde c, d ∈ R entao, pela Relacao de Euler, temos

eλt = e(c+di)t = ectedti = ect(cos(dt) + i sen(dt)).

Definicao 1.2. A funcao degrau unitario e definida por:

u(t) =

0 , se t ≤ 0

1 , se t > 0

Exemplo 1.1. A Transformada de Laplace da funcao degrau unitario e dada diretamente

da propriedade p3

L[u(t)] =1

s.

1.2 Fracoes Parciais

Quando estamos resolvendo uma EDO linear com coeficientes constantes, atraves da

transformada de Laplace, algumas vezes nos deparamos com uma equacao do tipo

L[y] =p(s)

q(s)


onde p(s) e q(s) sao polinomios e o grau de p(s) e menor que o de q(s).

Pelas propriedades da Transformada de Laplace conhecemos a solucao para

L[y] =u1

s− v1+

u2

s− v2+ . . .+

un

s− vn, (1.3)

onde uk, vk ∈ C para k = 1, 2, . . . , n, que e

y(t) = u1ev1t + u2e

v2t + . . .+ unevnt.

Considerando o caso onde as raızes de q(s) sao distintas podemos escrever L[y] na

forma da equacao 1.3, onde v1, v2, . . . , vn sao as raızes do polinomio q(s). Nesse caso

dizemos que a funcao L[y] foi decomposta em soma de fracoes parciais.

Exemplo 1.2. Decomponha F (t) em soma de fracoes parciais, onde

F (t) =10x+ 5

x2 − x− 6.

Sabendo que as raızes de x2−x−6 sao {−2, 3}, entao devemos determinar constantes

A1, A2 tais que possamos escrever F (t) na forma

F (t) =A1

x+ 2+

A2

x− 3

multiplicando ambos os membros por (x+ 2)(x− 3), temos

A1(x− 3) + A2(x+ 2) = 10x+ 5

logo, A1 + A2 = 10

−3A1 + 2A2 = 5

que possui solucao, A1 = 3 e A2 = 7, portanto

F (t) =3

x+ 2+

7

x− 3.


Proposicao 1.1. Se q(x) = sn + an−1sn−1 + . . . + a1s + a0, onde a0, a1, . . . , an−1 ∈ R

possui n raızes distintas e nao nulas {λ1, λ2, . . . , λn} entao a decomposicao de

L[y] =1

s q(x)

em soma de fracoes parciais, e

L[y] =A0

s+

A1

(s− λ1)+

A2

(s− λ2)+ . . .+

An

(s− λn)

donde A1, A2, . . . , An ∈ C sao dados por

A0 =1

(−λ1) . . . (−λn)

A1 =1

λ1(λ1 − λ2) . . . (λ1 − λn)...

An =1

λn(λn − λ1) . . . (λn − λn−1).

(1.4)

Demonstracao. Se q(x) = sn + an−1sn−1 + a1s+ a0 possui raızes λ1, λ2, . . . , λn, entao

q(x) = (s− λ1)(s− λ2) . . . (s− λn)

e para que L[y] seja decomposto em soma de fracoes parciais e necessario que exista

A0, A1, . . . , An ∈ C, tais que

L[y] =A0

s+

A1

(s− λ1)+

A2

(s− λ2)+ . . .+

An

(s− λn)

assim, multiplicando ambos os membros por s · q(x) temos

A0(s− λ1) . . . (s− λn) + A1s(s− λ2) . . . (s− λn) + . . .+ Ans(s− λ1) . . . (s− λn−1) = 1

assim, para s = 0, s = λ1, . . . , s = λn, temos

A0(−λ1)(−λ2) . . . (−λn) = 1

A1λ1(λ1 − λ2) . . . (λ1 − λn) = 1...

Anλn(λn − λ1) . . . (λn − λn−1) = 1

logo as equacoes dadas por sao validas.


Exemplo 1.3. Decomponha F (t) em soma de fracoes parciais, onde

F (t) =1

x3 − 3x2 − 4x.

Como x3− 3x2− 4x = x(x2− 3x− 4) e as raızes de x2− 3x− 4 sao {λ1 = 4, λ2 = −1},

entao devemos determinar constantes A0, A1 e A2 tais que possamos escrever F (t) na

forma

F (t) =A0

x+

A1

x− 4+

A2

x+ 1

pela proposicao 1.1 temos

A0 =1

(−λ1)(−λ2)=

1

(−4)(−(−1))= −1

4

A1 =1

λ1(λ1 − λ2)=

1

4(4 + 1)=

1

20

A2 =1

λ2(λ2 − λ1)=

1

(−1)(−1− 4)=

1

5

portanto

F (t) = −1

4

1

x+

1

20

1

x− 4+

1

5

1

x+ 1.

Notacao: Se A e um numero complexo, denotaremos por A o seu conjugado.

Algumas propriedades importantes dos numeros complexos

Se α, β ∈ C entao

1. α+ β = α+ β;

2. αβ = α β;

3. −α = −α;

4. α−1 = α−1.

Proposicao 1.2. Na proposicao 1.1, se λp = λp+1 entao Ap = Ap+1.


Demonstracao. Considere, sem perda de generalidade, que p = 1. Se λ1 = c + di = λ2,

onde c, d ∈ R, entao

(λ1 − λ2) = (λ2 − λ1)

e como cada um dos autovalores λk, (k > 2) da EDO ou e real ou e o conjugado de algum

outro autovalor da EDO (λk+1) e pode-se verificar que se λk e real, entao

(λ1 − λk) = (λ2 − λk)

e se λk nao e real, entao

(λ1 − λk)(λ1 − λk+1) = (λ1 − λk) (λ1 − λk+1)

= (λ1 − λk) (λ1 − λk+1)

= (λ2 − λk+1) (λ2 − λk)

= (λ2 − λk)(λ2 − λk+1)

assim, pela propriedades dos numeros complexos, temos A−11 = A1

−1= A−1

2 , pois

A−11 = λ1 (λ1 − λ2) (λ1 − λ3) . . . (λ1 − λn)

e

A−12 = λ2(λ2 − λ1)(λ2 − λ3) . . . (λ2 − λn)

portanto

A1 = A2.

Proposicao 1.3. Na proposicao 1.1, se λp e um numero real entao Ap tambem e um

numero real.

Demonstracao. Na proposicao 1.1, todas as raızes de q(x) = sn + an−1sn−1 + a1s + a0

sao distintas, logo se λj e uma raız de q(x) entao ou λj e um numero real ou λj e λj sao

raızes de q(x). Assim se λp e uma raiz real de q(x) entao dada uma raiz λj de q(x), j = p,

se λj ∈ R temos

(λp − λj) ∈ R


e se λj /∈ R entao λj tambem e uma raiz de q(x) assim,

(λp − λj)(λp − λj) ∈ R.

Logo

λp(λp − λ1) . . . (λp − λp−1)(λp − λp+1) . . . (λ1 − λn) ∈ R

e portanto

Ap =1

λp(λp − λ1) . . . (λp − λp−1)(λp − λp+1) . . . (λ1 − λn)∈ R.

1.3 Solucao e Analise da Estabilidade de EDO’s Li-

neares Com Coeficientes Constantes

Definicao 1.3. Sejam y(t) e r(t) duas funcoes de R em R, onde y(t) e n vezes dife-

renciavel, uma equacao diferencial ordinaria (EDO) e uma equacao na forma

f(y, y(1), . . . , y(n), t) = r(t).

Definicao 1.4. Uma Equacao Diferencial Ordinaria de ordem n e dita linear (EDOL)

quando pode ser escrita na forma

y(n) + pn−1(t)y(n−1) + . . .+ p1(t)y

(1) + p0(t)y = r(t). (1.5)

Se p1(t), p2(t), . . . , pn(t) forem funcoes constantes, a EDO e chamada linear com coe-

ficientes constantes reais. Se r(t) ≡ 0 dizemos que a EDO e homogenea e caso contrario

dizemos que ela e nao-homogenea.

Definicao 1.5. Dada a EDO linear com coeficientes constantes

y(n) + an−1y(n−1) + an−2y

(n−2) + ...+ a2y(2) + a1y

(1) + a0y = r(t)


chamamos de polinomio caracterıstico dessa EDO o polinomio

p(λ) = λn + an−1λn−1 + . . .+ a1λ+ a0 (1.6)

e chamamos de autovalores dessa EDO as raızes do seu polinomio caracterıstico.

Definicao 1.6. Definimos como Problema de Valor Inicial (PVI), a resolucao de uma

equacao diferencial de grau n, para t > 0, considerando as restricoes

y(k)(0) = y(k)0 , onde y

(k)0 ∈ R, para k = 0, 1, ..., n− 1

chamadas de condicoes iniciais.

Exemplo 1.4. Resolva a equacao diferencial ordinaria com coeficientes constantes

y′ − a0y = u(t)

com a condicao inicial y(0)=0, onde a0 = 0.

Aplicando a Transformada de Laplace em ambos os membros temos

L[y′]− a0L[y] = L[u(t)] ⇒ sL[y] + y(0)− a0L[y] = L[u(t)]

⇒ L[y](s− a0) = L[u(t)]

⇒ L[y] =1

s

1

s− a0

decompondo em soma de fracoes parciais temos

L[y] =1

a0

1

s− 1

a0

1

s− a0

logo, pela propriedade p3 da Transformada de Laplace apresentada nessa secao, a solucao

da EDO e

y(t) =1

a0− 1

a0ea0t.

Exemplo 1.5. Resolva a equacao diferencial ordinaria linear com coeficientes constantes

y′′ + 10y′ + 26y = 2u(t),

com condicoes iniciais y′(0) = y(0) = 0.


Aplicando, em ambos os membros a transformada de Laplace, e usando algumas pro-

priedades, temos

L[y′′] + 10L[y′] + 26L[y] = L[2u(t)] ⇒ s2L[y] + 10sL[y] + 26L[y] = 2L[u(t)]

⇒ L[y](s2 + 10s+ 26) =2

s

⇒ L[y] =2

s(s2 + 10s+ 26)

sabendo que as raızes de s2 + 10s+ 26 sao λ1 = −5 + i e λ2 = −5− i e decompondo L[y]

em uma soma de fracoes parciais temos

L[y] = 2

(1

26

1

s+−1 + 5i

52

1

s+ 5− i+−1− 5i

52

1

s+ 5 + i

)logo a solucao da EDO e

y(t) = 2

(1

26+−1 + 5i

52e(−5+i)t +

−1− 5i

52+ e(−5−i)t

)=

1

13+

(−1 + 5i

26

)e(−5+i)t +

(−1− 5i

26

)e(−5−i)t

portanto

y(t) =1

13+

1

13e−5t(cos(t)− 5sen(t)).

O grafico da solucao do exemplo 1.5 pode ser visto na figura 1.1.


0 0.5 1 1.5 2 2.5 30

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

Am

plitu

de

t (tempo)

Figura 1.1: Grafico da solucao do exemplo 1.5

Proposicao 1.4. Se a EDO linear com coeficientes constantes

y(n) + an−1y(n−1) + an−2y

(n−2) + ...+ a2y(2) + a1y

(1) + a0y = u(t) (1.7)

possui autovalores λ1, λ2, ..., λn (distintos e nao nulos), entao dadas as condicoes iniciais

y(0) = . . . = y(n−1)(0) = 0 a solucao geral da equacao e

y(t) = A0 + A1eλ1t + A2e

λ2t + ...+ Aneλnt (1.8)

onde A0 = 1/a0 e A1, A2, . . . , An sao numeros complexos tais que

A1 =1

λ1(λ1 − λ2)(λ1 − λ3) . . . (λ1 − λn)

A2 =1

λ2(λ2 − λ1)(λ2 − λ3) . . . (λ1 − λn)...

An =1

λn(λn − λ1)(λn − λ2) . . . (λn − λn−1).

Demonstracao. Aplicando a Transformada de Laplace a EDO dada pela equacao 1.7

temos

L[y(n)] + an−1L[y(n−1)] + an−2L[y

(n−2)] + ...+ a2L[y(2)] + a1L[y

(1)] + a0L[y] = L[u(t)]


donde, considerando as condicoes iniciais y(0) = . . . = y(n−1)(0) = 0, temos

L[y](sn + an−1sn−1 + . . .+ a2s

2 + a1s+ a0) = L[u(t)]

que implica que

L[y] =1

s(sn + an−1sn−1 + . . .+ a2s2 + a1s+ a0).

e se λ1, λ2, ..., λn (distintos e nao nulos) sao os autovalores da EDO, entao pela proposicao

1.1 temos que L[y] pode ser escrito na forma

L[y] =A0

s+

A1

(s− λ1)+

A2

(s− λ2)+ . . .+

An

(s− λn)

donde A1, A2, . . . , An sao numeros complexos tais que

A0 =1

(−λ1) . . . (−λn)

A1 =1

λ1(λ1 − λ2) . . . (λ1 − λn)...

An =1

λn(λn − λ1) . . . (λn − λn−1).

portanto

y(t) = A0 + A1eλ1t + A2e

λ2t + ...+ Aneλnt

e de

sn + an−1sn−1 + . . .+ a2s

2 + a1s+ a0 = (s− λ1)(s− λ2) . . . (s− λn)

temos que

a0 = (−λ1)(−λ2) . . . (−λn)

sendo assim

A0 =1

a0.

A partir da proposicao anterior e possıvel conhecer a solucao uma EDO com coefici-

entes constantes, a partir dos autovalores dessa EDO.


Corolario 1.1. Na proposicao 1.4, se λk = ck + dki e Ak = uk + vki, entao a equacao 1.8

pode ser escrita como

y(t) =1

a0+

n∑k=1

eckt(uk cos(dkt)− vk sen(dkt))

onde 1/a0 ∈ R e ck, dk, uk, vk ∈ R para todo k = 1, 2, . . . , n.

Corolario 1.2. Na proposicao 1.4, se a parte real de todos os autovalores da EDO dada

pela equacao 1.7 e nula e nenhum dos autovalores e nulo (λk = dki), entao sua solucao

pode ser escrita como

y(t) =1

a0+

n∑k=1

Ak cos(dkt)

onde 1/a0 ∈ R e Ak, dk ∈ R e dk = 0 para todo k = 1, 2, . . . , n.

Corolario 1.3. Na proposicao 1.4, se a parte imaginaria de todos os autovalores da EDO

dada pela equacao 1.7 e nula (λk = ck), entao sua solucao pode ser escrita como

y(t) =1

a0+

n∑k=1

Akeckt

onde 1/a0 ∈ R e Ak, ck ∈ R para todo k = 1, 2, . . . , n.

Exemplo 1.6. A figura 1.2 apresenta o grafico da solucao da EDO linear com coeficientes

constantes

y(4) + 20y(2) + 64y = u(t)

onde a parte real de todos os autovalores e nula. Nesse caso a funcao tem a forma dada

pelo corolario 1.2. Os autovalores sao {4i,−4i, 2i,−2i}.


0 2 4 6 8 100

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

Am

plitu

de

t (tempo)

Figura 1.2: EDO (exemplo 1.6) onde a parte real dos autovalores e nula

As proposicoes a seguir: 1.5 e 1.6 sao bem conhecidas. Suas demonstracoes podem ser

encontradas no apendice B.

Proposicao 1.5. Se a parte real de todos os autovalores da EDO 1.7 e negativa, entao

quando t→∞ sua solucao converge para 1/a0.

Exemplo 1.7. A figura 1.3 apresenta o grafico da solucao de duas EDOs lineares com

coeficientes constantes onde a parte real de todos os autovalores e negativa. A figura 1.3

(a) mostra o grafico da solucao da EDO

y(3) + 7y(2) + 14y(1) + 8y = u(t) (1.9)

que possui autovalores {−1,−2,−4}, com condicoes iniciais y0 = y′0 = y′′0 = 0. Enquanto

a figura 1.3 (b) mostra o grafico da solucao da EDO

y(5) + 8.2y(4) + 987.41y(3) + 5827.86y(2) + 66360.33y(1) + 234314.6y = u(t) (1.10)

que possui autovalores {−1 + 8i,−1− 8i,−1.1 + 30i,−1.1− 30i,−4}.


0 1 2 3 4 50

0.02

0.04

0.06

0.08

0.1

0.12

0.14

t

Am

plitu

de

(a) EDO da equacao 1.9

0 0.5 1 1.5 2 2.5 3 3.5 40

1

2

3

4

5

6x 10

−6

t

Am

plitu

de

(b) EDO da equacao 1.10

Figura 1.3: EDOs estaveis (a parte real de todos os autovalores e negativa)

Definicao 1.7. Uma EDO linear com coeficientes constantes e dita estavel quando a parte

real de todos os seus autovalores e negativa.

Exemplo 1.8. A EDO y′′ − 0.4y

′+ 1.04y = 1 e estavel.

Os autovalores dessa EDO sao λ1 = −0.2 + 1i e λ2 = −0.2− 1i que possuem a parte

real negativa, logo a EDO e estavel e converge quando t→∞ para

A0 =1

λ1λ2

=1

1.04= 0.9615384615.

Proposicao 1.6. Se a parte real de um dos autovalores da EDO 1.7 e positiva entao

quando t→∞ sua solucao diverge.

Exemplo 1.9. A figura 1.3 apresenta o grafico da solucao de duas EDOs lineares com

coeficientes constantes onde a parte real de um dos autovalores e positivo. A figura 1.3

(a) mostra o grafico da solucao da EDO

y(3) + 2y(2) − y(1) + 2y = u(t) (1.11)

que possui autovalores {1,−1,−2}, com condicoes iniciais y0 = y′0 = y′′0 = 0. Enquanto a

figura 1.3 (b) mostra o grafico da solucao da EDO

y(4) − 14y(3) + 50073y(2) − 380168y(1) + 400730144y = u(t) (1.12)

que possui autovalores {4− 100i, 3 + 200i, 3− 200i, 4 + 100i}.


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

3

t

Am

plitu

de

(a) EDO da equacao 1.11

0 0.2 0.4 0.6 0.8 1−1.5

−1

−0.5

0

0.5

1

1.5

2x 10

−7

t

Am

plitu

de

(b) EDO da equacao 1.12

Figura 1.4: EDOs instaveis (a parte real de um dos autovalores e positiva)

Definicao 1.8. Uma EDO linear com coeficientes constantes e dita instavel quando a

parte real de algum dos seus autovalores e negativa.

Exemplo 1.10. A EDO linear

y′′ + 0.4y′ + 1.04y = u(t)

sob condicoes y0 = y′0 = 0 e instavel.

Como essa EDO possui autovalores λ1 = 0.2 + 1i, λ2 = 0.2− 1i entao e instavel.

1.4 Pontos Crıticos da Solucao de EDO’s Lineares

com Coeficientes Constantes de Grau 2

Para uma funcao diferenciavel, ponto crıtico e um ponto onde a derivada sua derivada

e nula.

Proposicao 1.7. Dada a equacao diferencial linear com coeficientes constantes

y′′ + a1y′ + a0y = u(t)


onde a1, a0 ∈ R, possui autovalores λ1 = c + di e λ2 = c − di, onde c, d ∈ R e c < 0 e

d = 0 e u(t) = 0, se t ≤ 0 e u(t) = 1 se t > 0, com condicoes iniciais y(0) = 0 e y′(0) = 0.

Entao os pontos crıticos de y(t) sao

t =mπ

d, ∀ m ∈ N

e os valores de y(t) no pontos crıticos sao

y(t) =1

a0(1 + e

cmπd ), ∀ m ∈ N.

Demonstracao. Seja λ = c+ di, onde c = −a12e di =

√a21 − 4a0.

A partir das proposicoes 1.4 e 1.2 que a solucao para um problema desse tipo e da

forma

y(t) =1

a0+ veλt + veλt

onde v ∈ C. Denote k = λv, e pelas proepriedades dos numeros complexos, temos

y′(t) = keλt + keλt.

Sendo assim, como y′(0) = 0, temos k + k = 0 ⇒ k = qi, onde q ∈ R, temos

y(0) =1

a0+

qi

λ− qi

λ=

1

a0+

λqi

λλ− λqi

λλ=

1

a0+

λqi

a0− λqi

a0=

1

a0(1 + qi(λ− λ))

logo,

1 + qi[(c+ di)− (c− di)] = 1− 2dq = 0⇒ q =1

2d(1.13)

Assim,

y′(t) = 0 ⇒ keλt + keλt = 0

⇒ 1

2di eλt − 1

2di eλt = 0

⇒ 1

2di ect(cos(ct) + sen(dt))− 1

2di ect(cos(ct)− sen(dt)) = 0

⇒ 1

2di ect(2sen(dt)) = 0

⇒ sen(dt) = 0.


Portanto, y′(t) = 0 quando

t =mπ

d, ∀ m ∈ N.

Como k = λv, entao v =k

λ=

kλ

λλ=

kλ

a0e v =

kλ

a0, dessa forma podemos escrever

y(t) =1

a0+ veλt + veλt

=1

a0+ ect[(v + v)cos(dt) + (v − v)sen(dt)i]

=1

a0(1 + 2qect(d cos(dt) + c sen(dt))).

Sabendo que y′(mπd) = 0, para todo m ∈ N

y(mπ

d

)=

1

a0(1 + 2qd e

cmπd )

assim, pela equacao 1.13, os valores de y(t) no pontos crıticos sao

y(t) =1

a0(1 + e

cmπd ), ∀ m ∈ N.

Exemplo 1.11. Na figura 1.5 (a), temos os graficos de duas EDOs de grau 2

y′′ + 4y′ + 229y = 1 (1.14)

y′′ + 6y′ + 234y = 1 (1.15)

que possuem respectivamente autovalores {−2+15i,−2−15i}, {−3+15i,−3−15i} onde

as partes complexas dos seus autovalores sao iguais e sendo assim pela proposicao 1.7 o

tempos de maximo e mınimo de ambas as solucoes das equacoes sao iguais.

Na figura 1.5 (b) podemos ver os valores de maximo e mınimo da solucao das EDO’s

y′′ + 6y′ + 90y = 1 (1.16)

y′′ + 2y′ + 10y = 1 (1.17)

que possuem respectivamente autovalores {−3+9i,−3− 9i}, {−1+3i,−1− 3i}. A razao

entre a parte real e parte complexa dos autovalores das equacoes (−3/9 = −1/3) e a

mesma e pela proposicao 1.7 os valores de maximo e de mınimo da EDO devem ser os

mesmos.


0 0.5 1 1.5 20

1

2

3

4

5

6

7

8x 10

−3

exemplo 1

exemplo 2

(a) solucao da equacao 1.14 (exemplo 1) e da

equacao 1.15 (exemplo 2)

0 0.5 1 1.5 2 2.5 3 3.5 40

0.2

0.4

0.6

0.8

1

1.2

1.4

exemplo 1

exemplo 2

(b) solucao da equacao 1.16 (exemplo 1) e da

equacao 1.17 (exemplo 2)

Figura 1.5: Solucao das EDOs das equacoes 1.14, 1.14, 1.16 e 1.17

De acordo com a proposicao 1.7 o valor de maximo dos pontos crıticos da solucao de

uma EDO linear com coeficientes constantes de grau 2 com autovalores complexos (nao

reais) com a parte real negativa e determinado a partir da razao entre a parte real e a

parte complexa do autovalor da EDO, que pode ser visto tambem como a tangente do

angulo do autovalor.

Corolario 1.4. Se a equacao diferencial ordinaria linear com coeficientes constantes

y(n) + an−1y(n−1) + . . .+ a1y

(1) + a0y = 1

possui autovalores λ1, λ2, . . . , λn distintos e nao nulos e a equacao

x(n) + bn−1x(n−1) + . . .+ b1x

(1) + b0x = 1

possui autovalores pλ1, pλ2, . . . , pλn, onde p ∈ R, p = 0 entao para todo t ∈ R temos

x

(t

p

)=

1

pny(t).

Corolario 1.5. Dadas as condicoes da proposicao 1.4 (e pelo menos um dos autovalores

e complexo), se t0, t1, . . . sao os pontos de maximo e mınimo de y(t) entao t0/p, t1/p, . . .

sao os pontos de maximo e mınimo de x(t) e os valores de maximo e mınimo das equacoes

sao iguais a menos da constante pn.


1.5 Sistemas de Equacoes Lineares de Primeira Or-

dem

Uma EDO linear homogenea de ordem n com coeficientes constantes pode ser trans-

formada num sistema equivalente de equacoes de primeira ordem. E nesse caso resolver a

equacao caracterıstica da EDO sera equivalente a calcular os autovalores de uma matriz

especıfica.

Considere a seguinte EDO linear homogenea com coeficientes constantes de ordem n

y(n) + an−1y(n−1) + an−2y

(n−2) + . . .+ a2y(2) + a1y

(1) + a0y = 0

e defina x1 = y, x2 = y(1), . . . , xn−1 = y(n−2), xn = y(n−1), assim

x′1 = x2

x′2 = x3

...

x′n−1 = xn

e

x′n = −a0x1 − a1x2 − a2x3 − . . .− an−2xn−1 − an−1xn

dessa forma, se

A =

0 1 0 . . . 0

0 0 1 . . . 0...

......

. . ....

0 0 0 . . . 1

−a0 −a1 −a2 . . . −an−1

e x =

x1

x2

...

xn−1

xn

(1.18)

entao

x′ = Ax.


E evidente da definicao 1.5 que o polinomio caracterıstico da EDO e o polinomio

caracterıstico de A sao iguais, ou seja,

det(A− λI) = λn + an−1λn−1 + . . .+ λs+ a0. (1.19)

portanto os autovalores da EDO sao iguais aos autovalores da matriz A.

Logo, faz sentido definir as raızes do polinomio caracterıstico de uma equacao diferen-

cial ordinaria linear com coeficientes constantes como autovalores da EDO.

Definicao 1.9. Dada a equacao diferencial com coeficientes constantes

y(n) + an−1y(n−1) + an−2y

(n−2) + . . .+ a2y(2) + a1y

(1) + a0y = r(t)

dizemos que a matriz A definida em 1.18 e a sua matriz associada (ou companheira).

Exemplo 1.12. A EDO y′′+2y′+5y = 0 possui autovalores iguais aos da matriz D (sua

matriz associada), onde

D =

0 1

−5 −2

.

De fato,

D − λI =

−λ 1

−5 −2− λ

e

det(D − λI) = −λ(−2− λ) + 5 = λ2 + 2λ+ 5

portanto os autovalores de D sao iguais aos autovalores da EDO y′′ + 2y′ + 5y = 0, visto

que estes possuem o mesmo polinomio caracterıstico.

Se λ e um dos autovalores de uma EDO linear de coeficientes constantes reais, entao λ

tambem e autovalor dessa EDO, sendo assim os autovalores da matriz associada de uma

EDO possui apenas autovalores reais e/ou autovalores complexos conjugados.

Capıtulo 2

Metodo da Potencia

Resolver a equacao caracterıstica nao e a unica forma conhecida para calcular auto-

valores de uma matriz. Resolver uma equacao polinomial nao e facil para polinomios de

graus maiores que 4, e resolver uma equacao polinomial de alto grau atraves de calculo

numerico pode resultar em uma margem de erro maior do que se espera.

Nesse capıtulo sera apresentado o metodo da Potencia que permite calcular os auto-

valores que possui o maior valor absoluto de uma matriz. E atraves de uma adaptacao

calcular o autovalor que possui a maior parte real para assim conhecer a velocidade de

convergencia das EDOs lineares com coeficientes constantes estaveis.

Definicao 2.1. Dada uma matriz quadrada A de ordem n, dizemos que λ1 e um autovalor

dominante de A se

|λ1| ≥ |λ2| ≥ . . . ≥ |λn|.

Nas secoes seguintes serao apresentados dois algoritmos para calcular autovalores es-

pecıficos de uma matriz. O primeiro algoritmo (metodo da Potencia Classico) permite

calcular o autovalor dominante de uma matriz com um unico autovalor dominante e o

segundo algoritmo (metodo da Potencia Adaptado) permite calcular o par de autovalo-

res dominantes de uma matriz que possui somente um par de autovalores dominantes

25

Capıtulo 2. Metodo da Potencia 26

complexos conjugados. Como os autovalores de uma matriz associada a uma EDO sao

sempre reais ou pares de autovalores complexos conjugados, entao com esses dois algorit-

mos pode-se determinar os autovalores dominantes de qualquer matriz associada a uma

EDO linear.

Definicao 2.2. Dada um matriz A definimos exponencial dessa matriz como

eA =

(∞∑k=0

1

k!Ak

).

Proposicao 2.1. Se Av = λv entao Akv = λkv.

Demonstracao. Se Av = λv entao

A2v = A(Av) = Aλv = λAv = λ(λv) = λ2v

logo, pode-se mostrar por inducao que

Akv = λkv, ∀ k ∈ N.

Proposicao 2.2. Seja A uma matriz. Se λ e um autovalor de A, entao λ−1 e um autovalor

de A−1 (se A e inversıvel) e eλ e um autovalor de eA.

Demonstracao. Do teorema espectral temos que se λ e um autovalor da matriz A, entao

Av = λv

sendo assim, aplicando A−1 (considerando que A e inversıvel) em ambos os membros,

temos

A−1Av = A−1λv ⇒ v = λA−1v

⇒ λ−1v = A−1v.

Logo, λ−1 e autovalor de A−1.


Assim,

eAv =

(∞∑k=0

1

k!Ak

)v =

(∞∑k=0

1

k!Akv

)

=

(∞∑k=0

1

k!λkv

)=

(∞∑k=0

1

k!λk

)v = eλv.

Portanto, eλ e um autovalor de eA.

Corolario 2.1. Se λ e um autovalor da matriz A, entao eλi e um autovalor da matriz

eAi.

2.1 Metodo da Potencia para Matrizes com um Unico

Autovalor Dominante

Nessa secao sera feito o estudo do metodo da Potencia Classico (ver, por exemplo,

Watkins [13]) que trata do caso onde as matrizes possuem um unico autovalor dominante

, ou seja,

|λ1| > |λ2| ≥ . . . ≥ |λn|.

Proposicao 2.3. Seja A uma matriz de ordem n diagonalizavel e com um unico autovalor

dominante λ. Entao, para quase todo valor inicial x0 ∈ Rn dado, a sequencia de vetores

(xk) definida por

xk+1 =Axk

||Axk||, ∀ k = 0, 1, . . .

tende ao autovetor associado ao autovalor dominante de A.

Demonstracao. Considere que os autovalores estejam ordenados de modo que

|λ1| > |λ2| ≥ . . . ≥ |λn|

e sejam v1, v2, . . . , vn autovetores correspondentes, respectivamente, a λ1, λ2, . . . , λn. Agora

defina

x0 = c1v1 + c2v2 + . . .+ cnvn,


de modo que c1 = 0 (unica condicao sobre x0). Assim,

xk := Akx0 = Ak(c1v1 + c2 + v2 + . . .+ cnvn)

= Akc1v1 + Akc2v2 + . . .+ Akcnvn

e pela proposicao 2.1 temos

xk = Akx0 = λk1c1v1 + λk

2c2v2 + . . .+ λkncnvn

logo

xk = λk1(c1v1 +

(λ2

λ1

)k

c2v2 + . . .+

(λn

λ1

)k

cnvn)

desde que λ1 = 0.

Assim, se λ1 e o unico autovalor dominante de A, entao |λm/λ1| < 1, o que implica

que

(λm/λ1)k → 0 ∀ m = 2, 3, . . . , n

quando k →∞.

Agora defina

xk =xk

||xk||.

Sendo assim, quando k →∞ temos

xk =λk1c1v1

||λk1c1||||v1||

= ± v1||v1||

portanto xk = Akx0

||Akx0|| tende ao autovetor de A associado a autovetor que possui o maior

valor absoluto.

Se A e uma matriz quadrada de ordem n diagonalizavel e com um unico autovalor

dominante λ, entao e possıvel determinar λ pelo algoritmo 2.1.


Algoritmo 2.1: Metodo da Potencia (Classico)

Entrada: A, x0, numInt (numero de iteracoes), tol (tolerancia).

Saıda: λ (Autovalor Dominante de A).

inıcio

x0 ← x0/|x0|;

λ0 ← xt0Ax0;

para k=1 ate numInt faca

xk ← Axk−1;

xk ← xk/|xk|;

λk ← xtkAxk;

se |λk − λk−1| < tol entao

break; /* Sair do Laco */

fim se

fim para

λ← λk;

fim

De acordo com a proposicao 2.3, no algoritmo 2.1 o vetor xk final representa uma apro-

ximacao para o autovetor associado ao autovalor dominante de A. Sendo assim escolhendo

o numero de iteracoes (numInt) e a tolerancia (tol) adequados podemos escrever

Axk = λxk

donde

λ =xtkAxk

xtkxk

.

E muito raro que o valor de c1 da proposicao 2.3 seja nulo, se o vetor x0 no algoritmo 2.1

for escolhido de forma aleatoria e ainda nesse caso, gracas a erros do calculo computacional

o metodo ainda funciona, mas requer um aumento no numero de iteracoes do metodo.


Exemplo 2.1. Considere a matriz

A =

0 2

10 −8

,

que tem como autovalores e respectivos autovetores

λ1 = −10, v1 =

1/√26

−5/√26

≈ −0.196116135138184

0.980580675690920

e

λ2 = 2, v2 =

1/√2

1/√2

≈ 0.707106781186547

0.707106781186547

donde λ1 e o autovalor dominante.

O metodo da Potencia aplicado a matriz A com numInt = 9 e tol = 10−5 e

x0 = [0.657818756480074, 0.753176263316237]t

escolhido aleatoriamente, retorna o valor aproximado para o autovalor dominante da

matriz A igual a

λ = −10.000040068812652.

Os valores de xk e λk a cada iteracao k que serao retornados pelo metodo da Potencia,

dado pelo algoritmo 2.1, com esse valores iniciais podem ser vistos na tabela 2.1.

De acordo com a tabela 2.1, podemos observar a rapida convergencia para o valor

exato −10. A velocidade da convergencia esta associada ao valor absoluto da razao entre

os autovalores ∣∣∣∣λ2

λ1

∣∣∣∣ = 2

10= 0, 2

quanto mais proximo de 1 for essa razao mais lenta sera a convergencia do metodo. Isso

ocorre pois na proposicao 2.3 e necessario que(λi

λ1

)k

→ 0, quando k → +∞

para todo i = 2, 3, . . . , n, e essa convergencia e mais rapida quando o valor absoluto da

razao entre os autovalores e o autovalor dominante for mais proximo de 0.


k xk λk

0 [ 0.657818756480074, 0.753176263316237]t 1.407245806355873

1 [ 0.938786002137701, 0.344500859491398]t 2.931504277856682

2 [ 0.103336593054350, 0.994646443986968]t -6.681171890504509

3 [ 0.276140582947024,−0.961117255307425]t -10.574812777529857

4 [−0.180904934123345, 0.983500587091756]t -9.873228545495890

5 [ 0.199191173521070,−0.979960650430056]t -10.024977154694001

6 [−0.195502433185524, 0.980703216380236]t -9.994988769370785

7 [ 0.196238927810593,−0.980556109160383]t -10.001001619931117

8 [−0.196091578694559, 0.980585586659867]t -9.999799650919769

9 [ 0.196121046510547,−0.980579693403656]t -10.000040068812652

Tabela 2.1: Iteracoes do metodo da Potencia aplicado a matriz do exemplo 2.1

2.2 Metodo da Potencia para Matrizes com um Par

de Autovalores Dominantes Complexos Conjuga-

dos

A matriz associada de uma EDO linear com coeficientes constantes possui apenas

autovalores reais ou complexos conjugados. Autovalores complexos conjugados possuem

o mesmo valor absoluto, sendo assim e comum que a matriz associada de uma EDO linear

com coeficientes constantes possua um par de autovalores complexos conjugados.

O metodo da Potencia Classico nao permite calcular esses autovalores, visto que a

convergencia do metodo acontece apenas quando a matriz possui um unico autovalor do-

minante. Esse e um problema bem incomum e levou ao desenvolvimento de um algoritmo

especıfico, visto que o metodo da Potencia classico nao funciona.

Assim, diferentemente das secoes anteriores, nessa secao serao apresentados resultados

novos (pelo menos que nao foram encontrados na literatura) como um algoritmo para cal-


cular o autovalor dominante de uma matriz que possui um par de autovalores dominantes

e estes sao complexos conjugados.

Proposicao 2.4. Sejam A uma matriz real diagonalizavel com um par de autovalores

complexos conjugados dominantes, λ = c+di e λ, onde c, d ∈ R, e respectivos autovetores

associados v = v1 + v2i e v, onde v1, v2 ∈ Rn. Se

u0 = g1v1 + g2v2

onde g1, g2 ∈ R∗, entao existem α, β ∈ R definidos a partir de

A2u0 = αu0 + βAu0

tais que,

λ =β +

√β2 + 4α

2.

Demonstracao. Seja λ um autovalor dominante de A, donde λ = c+ di, com c, d ∈ R e

v = v1 + v2i o autovetor associado a λ, com v1, v2 ∈ Rn. Dessa forma Av = λv, logo

A(v1 + v2i) = (c+ di)(v1 + v2i)

= (cv1 − dv2) + (cv2 + dv1)i

que pode ser escrito como

A[v1 v2

]=[v1 v2

] c d

−d c

.

Agora defina u1 = Au0 e u2 = A2u0, assim temos

u1 = Au0 = A(g1v1 + g2v2) = g1Av1 + g2Av2 = g1(cv1 − dv2) + g2(cv2 + dv1)

ou seja,

u1 = (g1c+ g2d)v1 + (g2c− g1d)v2

denotando g3 = g1c+ g2d e g4 = g2c− g1d temos, pelo mesmo raciocınio, que

u2 = (g3c+ g4d)v1 + (g4c− g3d)v2.


logo u0, u1, u2 ∈ span{v1, v2}. E definindo G como

G =

g1 g1c+ g2d

g2 (g2c− g1d)

temos [

u0 u1

]=[v1 v2

]G.

Caso λ ∈ R, teremos v2 = 0, logo u0 = g1v1 e u1 = (g1c)v1 e portanto

span{u0, u1} = span{u0} = span{v1} = span{v1, v2}

e caso λ /∈ R pode-se verificar que a matriz G e inversıvel e portanto tambem teremos

span{u0, u1} = span{v1, v2}.

De todo modo teremos u2 ∈ span{u0, u1}. Portanto existem constantes α e β, que

satisfazem

u2 = αu0 + βu1 (2.1)

ou seja,

A2u0 = αu0 + βAu0.

Assim, escreva

A[u0 u1

]=

[Au0 Au1

]=

[u1 u2

]=

[u1 αu0 + βu1

]ou seja,

A[u0 u1

]=

[u0 u1

] 0 α

1 β

e com isso temos[

v1 v2

]G

0 α

1 β

= A[v1 v2

]G =

[v1 v2

] c d

−d c

G


logo c d

−d c

G = G

0 α

1 β

assim as matrizes c d

−d c

e

0 α

1 β

tem os mesmos autovalores, visto que G e inversıvel, e calculando diretamente seus auto-

valores, temos que, respectivamente, seus autovalores sao

c± di eβ ±

√β2 + 4α

2.

Portanto, temos que

λ =β +

√β2 + 4α

2.

Alguns resultados usados na demonstracao da proposicao anterior podem ser encon-

trados em Horn e Johnson [6].

Proposicao 2.5. Seja A uma matriz real diagonalizavel com um par de autovalores com-

plexos conjugados dominantes, λ /∈ R e λ cujo autovetores associados sao, respectiva-

mente, v = v1 + v2i e v com v1, v2 ∈ Rn. Entao para quase todo valor inicial x0 ∈ Rn

dado, a sequencia (xk) definida por

xk+1 =Axk

||Axk||

satisfaz

infαk,βk∈R

{||xk − αkv1 − βkv2||} = 0.

Demonstracao. A demonstracao dessa proposicao e semelhante a demonstracao da pro-

posicao 2.3.


A condicao, quase todo, em x0 na proposicao 2.5 e equivalente ao da proposicao 2.3.

As proposicoes 2.4 e 2.5 (juntas) permitem desenvolver um algoritmo para calcular

os autovalores dominantes de uma matriz que possui um par de autovalores complexos

conjugados dominantes. Com esse algoritmo sera possıvel resolver o problema de calcular

os autovalor dominante da matriz associada a uma EDO. Como a proposicao 2.5 e uma

leve modificacao da proposicao 2.3 (metodo da Potencia) daremos a esse algoritmo o nome

de metodo da Potencia Adaptado.


Algoritmo 2.2: Metodo da Potencia (Adaptado)

Entrada: A, x0, numInt, tol.

Saıda: λ (Autovalor Dominante de A).

inıcio

N0 ← 1/|x0|;

x0 ← x0N0;

x1 ← Ax0;

N1 ← 1/|x1|;

x1 ← x1N1;

λ0, λ1, e0, e1 ← 0;

para k=2 ate numInt faca

xk ← Axk−1;

Nk ← 1/|xk|;

xk ← xkNk;

v0 ← uk−2;

v1 ← uk−1Nk−2;

v2 ← ukNk−1Nk−2;

resolva: v2 = a v0 + b v1;

λk ← (b+√b2 + 4a)/2;

ek ← |λk − λk−1|;

se ek < tol entao

break; /* Sair do Laco */

fim se

fim para

λ← λk;

fim


Exemplo 2.2. Considere a matriz

A =

0 1 0

0 0 1

−60 −32 −7

,

que tem como autovalores e respectivos autovetores

λ1 = λ2 = −2 + 4i, v1 = v2 =

−0.0292 + 0.0390i

−0.0975− 0.1949i

0.9747

, λ3 = −3, v3 =

−0.1048

0.3145

−0.9435

donde λ1 e λ2 sao os unicos autovalores dominantes de A e sao complexos conjugados.

O metodo da Potencia aplicado a matriz A com numInt = 19, tol = 10−4 e

x0 = [0.8742, 0.4856, 0.0043]t

escolhido aleatoriamente, retorna o valor aproximado para o autovalor dominante da

matriz A igual a

λ = −2.0001 + 4.0021i.

Os valores de xk e λk a cada iteracao k que serao retornados pelo metodo da Potencia,

dado pelo algortimo 2.2, com esses valores iniciais podem ser vistos na tabela 2.2.

E possıvel diminuir o tempo computacional desse metodo, deixando de calcular o

valor de λk em cada passo, mas nesse caso teremos apenas um criterio de parada para o

algoritmo.


k xk λk ek

1 [0.8742, 0.4856, 0.0043]t 0 0

2 [0.0071, 0.0001,−1.0000]t 0 0

3 [0.0000,−0.1505, 0.9886]t -3.2857 + 4.2878i 5.4019

4 [−0.0646, 0.4241,−0.9033]t -2.4719 + 3.5182i 1.1201

5 [0.1205,−0.2566,−0.9590]t -1.8773 + 3.5258i 0.5946...

......

...

15 [0.0068,−0.2733, 0.9619]t -2.0099 + 3.9964i 0.0145

16 [−0.1447, 0.5090, 0.8485]t -2.0008 + 3.9930i 0.0097

17 [0.0375, 0.0625,−0.9973]t -1.9960 + 3.9974i 0.0065

18 [0.0215,−0.3426, 0.9393]t -1.9972 + 4.0016i 0.0044

19 [−0.1052, 0.2884, 0.9517]t -2.0001 + 4.0021i 0.0029

Tabela 2.2: Iteracoes do metodo da Potencia Adaptado aplicado a matriz do exemplo 2.2

2.3 Metodo da Potencia para o Calculo do Autovalor

com Maior Valor Absoluto de uma Matriz

Proposicao 2.6. Seja A uma matriz quadrada de ordem n. O metodo da Potencia

aplicado a eA converge para eλ, onde λ e o autovalor de A que possui a maior parte real.

Demonstracao. Sejam λ1 = c1 + d1i, λ2 = c2 + d2i, . . . , λn = cn + dni os autovalores de

A, onde ck, dk ∈ R ∀k = 1, 2, ..., n e µ1, µ2, . . . , µn os autovalores de eA, assim o metodo

da Potencia aplicado a eA converge para µ1 se

|µ1| > |µ2| ≥ |µ3| ≥ . . . ≥ |µn|

mas, pela proposicao 2.2, os autovalores de eA sao

µ1 = eλ1 , µ2 = eλ2 , . . . , µn = eλn


sendo assim o metodo da Potencia aplicado a eA converge para µ1 se

|eλ1| > |eλ2 | ≥ |eλ3 | ≥ . . . ≥ |eλn|

mas

|eλk | = |eck+dki| = |eck(cos(dk) + sen(dk)i)| = |eck | |cos(dk) + sen(dk)i| = |eck |

logo, o metodo da Potencia aplicado a eA converge para µ1 = eλ1 se

|ec1 | > |ec2 | ≥ |ec3 | ≥ . . . ≥ |ecn|

ou seja, se

c1 > c2 ≥ c3 ≥ . . . ≥ cn.

Portanto o metodo da Potencia aplicado a eA converge para µ1 = eλ1 , onde a parte

real de λ1 e a maior dentre os autovalores de A, ou seja, o metodo da Potencia aplicado

a eA converge para o autovalor de eA que esta associado ao autovalor de A que possui a

maior parte real.

Corolario 2.2. Se A possui apenas autovalores reais ou complexos conjugados entao eA

tambem possui apenas autovalores reais ou complexos conjugados.

Exemplo 2.3. A matriz

A =

0 1 0 0

0 0 1 0

0 0 0 1

−200 −260 −69 −10

,

possui autovalores λ1 = −1, λ2 = λ3 = −2 + 6i e λ4 = −5.

Ao aplicarmos o metodo da Potencia (Adaptado) a matriz A com

numInt = 17, tol = 10−8 e x0 = [0.1808, 0.2456, 0.7964, 0.8558]t


obteremos o seguinte valor de λ

λ = −2.001517020121187 + 6.105585388659423i

e λ2 = λ3 = −2+6i sao os unicos autovalores dominantes deA e sao complexos conjugados.

Mas caso o interesse seja determinar o autovalor de A que possui a maior parte real, pode-

se aplicar o metodo da Potencia (Adaptado) a eA.

Aplicando o metodo da Potencia (Adaptado) a eA com os mesmos valores de entrada

teremos como retorno

µ = 0.367877861118490,

donde µ = eλ, logo

λ = ln(µ) = −1.000004295038453

que e uma boa aproximacao, visto que o autovalor de A que possui a maior parte real e

λ1 = −1.

2.4 Tempo de Acomodacao de EDO’s Lineares Com

Coeficientes Constantes

Dada uma EDO linear com coeficientes constantes estavel sem autovalores repetidos

(no maximo complexos conjugados) conhecemos a sua solucao que de acordo com o Co-

rolario 1.1

y(t) = A0 +n∑

k=1

eckt(ukcos(dkt) + vksen(dkt)).

Como a EDO e estavel sabemos que os valores c1, c2, . . . , cn sao todos negativos, e que

y(t) converge para A0. A medida que o valor de t aumenta o autovalor que possui a maior

parte real da EDO tera mais influencia sobre o valor de y(t), visto que ecit, converge mais

rapidamente que ecjt, quando ci < cj < 0.

Sendo assim, se as partes reais dos autovalores estao ordenadas de modo que

c1 > c2 ≥ . . . ≥ cn


podemos considerar a funcao f1(t) = A0+ec1t(u1cos(d1t)+v1sen(d1t)) como uma simples

aproximacao para y(t). A funcao f1(t) e uma aproximacao para y(t) a medida que o valor

de t aumenta, mas nao e tao simples determinar o valor de u1cos(d1t)+v1sen(d1t). Sendo

assim considere a funcao fa(t) = A0−A0ec1t para estudar a estabilidade assintotica de y(t)

visto que nao precisamos calcular o valor de u1cos(d1t) + v1sen(d1t) que e determinado

em funcao dos autovalores, alem de satisfazer y(0) = fa(0) = 0. Alem disso essa funcao

tambem e solucao do PVI

y′ − a0y = u(t), y0 = 0,

veja o exemplo 1.4.

Definicao 2.3. Dada uma equacao diferencial ordinaria linear com coeficientes constantes

estavel com autovalores distintos

y(n) + an−1y(n−1) + . . .+ a1y

(1) + a0y = u(t)

a funcao de acomodacao de y e a funcao

fa(t) =1

a0(1− ect)

onde c e a parte real do autovalor que possui a maior parte real dentre os autovalores da

EDO.

A funcao de acomodacao de um EDO e uma das aproximacoes mais simples e de facil

manejo que podemos ter da solucao exata da EDO, alem de ser uma funcao monotona e

convergente. Na engenharia a funcao de acomodacao tambem e conhecida como funcao

de estabelecimento ou de assentamento.

Pode-se observar tambem que a solucao e a funcao de acomodacao da EDO

y′ − a0y = u(t) sao iguais.

Na tabela 2.4 temos quatro EDO’s lineares com coeficientes constantes (sendo um de

grau 2, uma de grau 3 e duas de grau 4) e seus respectivos autovalores.


Ex EDO Autovalores

1 y′′ + 6y′ + 58y = 1 -3 + 7i, -3 - 7i

2 y′′′ + 8y′′ + 45y′ + 116y = 1 -2 + 5i, -2 - 5i -4

3 y′′′′ + 24y′′′ + 875y′′ + 7594y′ + 41114y = 1 -5+6i, -5-6i, -7+25i, -7-25i

4 y′′′′ + 23y′′′ + 286y′′ + 2140y′ + 5800y = 1 -4+10i,-4-10i, -5, -10

Tabela 2.3: EDO’s lineares com coeficientes constantes para estudar o tempo de aco-

modacao

Na figura 2.1 temos graficos com a solucao das quatro equacoes diferenciais da tabela

2.4 e suas respectivas funcoes de acomodacao.

Definicao 2.4. O tempo de acomodacao com precisao p, com p ∈ (0, 1), denotado por ta,

da EDO linear com coeficientes constantes estavel

y(n) + an−1y(n−1) + . . .+ a1y

(1) + a0y = u(t)

e o valor de t que satisfaz

fa(t) = p1

a0

onde fa e a funcao de acomodacao da EDO.

Assim como a funcao de acomodacao, o tempo de acomodacao tambem e conhecido

como tempo de estabelecimento ou de assentamento. Pela definicao de tempo de aco-

modacao, se estivermos interessados em determinar o tempo de acomodacao de uma EDO

com precisao de 99%, precisamos determinar o valor de t, tal que a funcao de acomodacao

nesse valor e 99% do valor limite, 1/a0.

Proposicao 2.7. O tempo de acomodacao da EDO estavel

y(n) + an−1y(n−1) + . . .+ a1y

(1) + a0y = u(t)

e

ta =ln(1− p)

c(2.2)

onde c e a parte real do autovalor que possui a maior parte real dentre os autovalores da

EDO.


0 0.5 1 1.5 2 2.5 30

0.005

0.01

0.015

0.02

0.025

(a) exemplo 1

0 1 2 3 4 50

0.001

0.002

0.003

0.004

0.005

0.006

0.007

0.008

0.009

0.01

(b) exemplo 2

0 0.5 1 1.5 20

0.5

1

1.5

2

2.5

3x 10

−5

(c) exemplo 3

0 0.5 1 1.5 20

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8x 10

−4

(d) exemplo 4

Figura 2.1: Analise do tempo de acomodacao das EDO’s da tabela 2.4

Demonstracao. Pela definicao de tempo de acomodacao temos que

fa(ta) = p1

a0

sendo c a parte real do autovalor que possui a maior parte real dentre os autovalores da

EDO, temos

fa(ta) =p

a0⇒ 1

a0− 1

a0ecta =

p

a0⇒ 1− ecta = p

⇒ ecta = 1− p

⇒ cta = ln(1− p)


logo,

ta =ln(1− p)

c.

Iremos chamar de Regime Transitorio da EDO lineares com coeficientes constantes

estaveis o conjunto de valores no tempo que sao menores que tempo de acomodacao e de

Regime Permanente ou Estacionario o conjunto de valores no tempo que sao maiores que

o tempo de acomodacao.

Capıtulo 3

Metodos Splines para o Calculo

Aproximado da Solucao de EDOs

Nesse capıtulo sera descrito Spline polinomial (Burden e Faires [3]) e serao apresen-

tados metodos numericos para o calculo aproximado da solucao de EDOs baseado em

tecnicas de interpolacao por splines que chamaremos de metodo da Sub-particao, metodo

da Suavidade e metodo da Suavidade Adaptado.

Esses metodos sao capazes de aproximar a solucao de EDO gerais (incluindo equacoes

nao lineares e nao estavel). Porem mesmo no caso linear estavel eles parecem ser competi-

tivos com os metodos normalmente usados, como o calculo da solucao exata e os metodos

de Runge-Kutta.

Quando a EDO a ser estudada for linear com coeficientes constantes e estavel esses

metodos determinarao um spline de aproximacao para a regiao da solucao da EDO menos

estavel, conhecido como regime transitorio, visto que no regime permanente a solucao

pode ser aproximado por uma funcao constante.

Por fim, serao apresentados resultados numericos comparativos com o metodo de Euler

(Boyce e DiPrima [2], Mathews [9]), metodo numerico classico para resolucao numerica

45

Capıtulo 3. Metodos Splines para o Calculo Aproximado da Solucao de EDOs 46

de EDO’s, em tempo computacional e erro relativo.

Com excecao dessa primeira secao, que trata sobre Spline de modo geral, todas as

outras secoes apresentam resultados novos (que nao foram encontrados na literatura), ou

seja, o metodo da Sub-particao, o metodo da Suavidade, o metodo da Suavidade Adaptado

que permitem o calculo aproximado da solucao de EDO’s lineares sao apresentados pela

primeira vez nesse trabalho.

3.1 Spline

A tecnica de interpolacao por Spline e usada quando se deseja passar por um conjunto

de pontos dados no plano uma funcao contınua e suave (sem bicos). Uma discussao

completa sobre esse tema pode ser encontrado em Ahlberg, Nilson, and Walsh [1].

Definicao 3.1. Sejam xj, yj ∈ R para j = 1, 2, . . . , nk com n, k ∈ N∗ e x1 > x2 > . . . >

xnk. Um Spline de grau n e ordem k e uma funcao polinomial por partes de classe Ck

que satisfaz

S(xj) = yj ∀j = 1, 2, . . . , nk.

O problema geral que leva ao uso de spline e interpolar por um conjunto de pontos

dados no plano uma funcao que e polinomial por partes.

Uma forma de determinar os coeficientes de cada um dos polinomios

Si(x), i = 1, 2, . . . , n− 1

e montar um sistema linear baseado nas condicoes que definem o spline S(x). Considere

uma situacao onde se deseja determinar um spline de grau 3 e ordem 1 que interpola os

pontos no plano, p1 = (x1, y1), p2 = (x2, y2), . . . , pn = (xn, yn). Para isso defina S0(x) = 0

e

Sk(x) = akx3 + bkx

2 + ckx+ dk, para k = 1, 2, . . . , n− 1.


Em cada intervalo [xk, xk+1] vamos determinar constantes ak, bk, ck, dk tais que possamos

definir uma funcao que passe pelos pontos pk = (xk, yk) dados, e seja de classe C2.

Essas quatro variaveis, para cada k = 1, 2, . . . , n−1 sao determinadas por um sistema

de equacao

Sk(xk) = yk

S ′k(xk) = S ′

k−1(xk)

S ′′k (xk) = S ′′

k−1(xk)

Sk(xk+1) = yk+1

(3.1)

sendo assim pode-se definir a funcao (spline), S(x), que passa pelos pontos dados

S(x) = Si(x), se xi < x ≤ xi+1.

Nesse trabalho iremos chamar as tres primeiras linhas do sistema 3.1 de condicoes de

suavidade e a ultima linha de condicao de interpolacao (a primeira linha tambem pode

ser considerada uma condicao de interpolacao).

Exemplo 3.1. Dados os pontos p1 = (1, 3), p2 = (2, 5), p3 = (4, 2) e p4 = (6, 3) determi-

nar uma funcao que passa por esses pontos e seja de classe C2 no intervalo [1, 6].

Substituindos os pontos dados no sistema de equacoes 3.1 para k = 1, 2, 3, temos tres

sistemas de equacoes lineares e resolvendo cada um temos

S(x) =

2x3 − 6x2 + 6x + 1, se 1 ≤ x ≤ 2

−4.875x3 + 35.25x2 − 76.5x + 56, se 2 < x < 4

18.875x3 − 249.75x2 + 1063.5x − 1464, se 4 ≤ x ≤ 6

A determinacao e solucao desses sistemas foram calculadas numericamente o software

Matlab.

O grafico de S(x) pode ser visto na figura 3.1.


1 2 3 4 5 6−35

−30

−25

−20

−15

−10

−5

0

5

10

15S

(x)

x

(a) Grafico de S(x)

1 2 3 4 5 6

2

2.5

3

3.5

4

4.5

5

S(x

)

x

(b) Grafico de S(x) ampliado

Figura 3.1: Grafico da solucao do exemplo 3.1

3.2 Metodo da Sub-particao

Nessa secao sera apresentado o metodo da Sub-particao que permite determinar um

spline de grau m e classe C1, como solucao aproximada de uma Equacao Diferencial

Ordinaria. Esse metodo foi criado especificamente para esse trabalho.

Considere a seguinte EDO geral de ordem n

f(y, y(1), . . . , y(n), t) = r(t)

estamos interessados em determinar uma funcao (spline) que seja uma aproximacao para

solucao exata dessa EDO, y(t). O processo se inicia definindo o intervalo [0, xf ] onde

sera feita aproximacao quando a EDO que se deseja aproximar a solucao for linear com

coeficientes constantes pode-se tomar xf = ta. O intervalo [0, xf ] deve ser dividido, em s

intervalos, onde s ∈ N, de forma que o comprimento de cada um deles seja o mesmo (por

questao de simplicidade). Definindo h como sendo o comprimento de cada intervalo temos

sh = xf . Sendo assim, temos s intevalos: I1 = [0, h], I2 = [h, 2h], . . . , Is = [(s − 1)h, sh]

e em cada um desses intervalos o spline de aproximacao S(t) sera um polinomio de grau

m, onde m > 2.

Para determinar os polinomios Sk(x), para k = 1, 2, . . . , s, que definem o spline S(t)


resolve-se os seguintes sistemas para k = 1, 2, . . . , s

Sk(xk) = Sk−1(xk)

S ′k(xk) = S ′

k−1(xk)

f(Sk(xk + ( 1m−1

)h), . . . , S(n)k (xk + ( 1

m−1)h), xk + ( 1

m−1)h) = r(xk + ( 1

m−1)h)

f(Sk(xk + ( 2m−1

)h), . . . , S(n)k (xk + ( 2

m−1)h), xk + ( 2

m−1)h) = r(xk + ( 2

m−1)h)

......

...

f(Sk(xk + (m−2m−1

)h), . . . , S(n)k (xk + (m−2

m−1)h), xk + (m−2

m−1)h) = r(xk + (m−2

m−1)h)

f(Sk(xk+1), . . . , S(n)k (xk+1), xk+1) = r(xk+1)

(3.2)

(Sistema de equacoes para o metodo da Sub-particao)

donde xk = (k − 1)h e S0(t) satisfaz as condicoes iniciais, S0(0) = y0 e S ′0(0) = y′0.

Para concluir define-se o spline

S(t) = Sk(t), se xk−1 < t ≤ xk.

A primeira e a segunda linha do sistema de equacoes 3.2 exige do spline que ele seja

contınuo e que sua derivada seja contınua no intervalo [0, ta] e a ultima linha do sistema

3.2 exige que o spline satisfaca a EDO no pontos xk para k = 1, 2, . . . , s. As outras linhas

exigem que o spline satisfaca a EDO em m− 2 outros pontos para cada k = 1, 2, . . . , s.

O nome de Sub-particao foi dado, por mim, para esse metodo pelo fato desse metodo

exigir pouco do spline quanto a sua suavidade (apenas as duas primeiras linhas do sistema

3.2), e exigir que o spline satisfaca a EDO em muitos pontos, visto que cada intervalo

[xk, xk+1], foi dividido na maior quantidade de partes possıvel, a depender apenas do grau

do spline no intervalo. Nas duas primeiras linhas do sistema 3.2 sempre temos equacoes

lineares. Ja nas outras linhas teremos equacoes lineares se a EDO f(y, y(1), . . . , y(n), t) =

r(t) for linear e equacoes nao lineares caso contrario.

Como o metodo da Sub-particao usa apenas duas condicoes iniciais ele nao e ideal

para determinar uma funcao de aproximacao, porem e boa para conhecer a acomodacao

da solucao das EDO lineares com coeficientes constantes.


Algoritmo do metodo da Sub-particao

1. Escolha um valor para xf , caso a EDO em questao for linear com coeficientes cons-

tantes pode-se determinar que xf = ta, e calcular ta (tempo de acomodacao) usando

o Metodo da Potencia para calcular a parte real do autovalor que possui a maior

parte real da EDO.

2. Escolha o valor de s de forma conveniente e calcule o valor de h pela equacao

sh = xf .

3. Defina xk = (k − 1)h.

4. Defina o polinomio S0(x) de maneira que este satisfaca as condicoes iniciais para a

primeira e segunda derivada.

5. Determine os coeficientes dos polinomios Sk(x), para k = 1, 2, . . . , s atraves do

sistema de equacoes 3.2.

6. Defina o spline S(t) = Sk(t) se xk−1 < t ≤ xk.

No restante dessa secao tratatemos apenas de EDO lineares com coeficientes constan-

tes.

Exemplo 3.2. Calcule a solucao aproximada usando o metodo da Sub-particao do se-

guinte problema de valor inicial

y′′ + 4y′ + 13y = 1 (3.3)

sujeito a y(0) = 0 e y′(0) = 0.

Os autovalores desse problema, λ1 = −2 + 3i e λ2 = −2− 3i, podem ser rapidamente

calculados.

Como se trata de uma EDO linear com coeficientes constantes entao vamos considerar

xf = ta. Tomando s = 3 (tres polinomios), e escolhendo o tempo de acomodacao com


precisao de 99% podemos determinar o valor de h, assim que for calculado o valor de ta

pela equacao 2.2, ou seja,

ta = −1

2ln(0, 01) ≈ 2.3026 ⇒ h =

2.3025

3≈ 0.7675.

Agora defina

Sk(t) = akt3 + bkt

2 + ckt+ dk, k = 1, 2, 3

e S(t), a aproximacao (spline) da solucao da EDO, como sendo

S(t) =

0 se t ≤ 0

S1(t) se 0 < t ≤ h

S2(t) se h < t ≤ 2h

S3(t) se 2h < t ≤ 3h

1/13 se t > 3h

Para determinar os parametros a1, b1, c1, d1 e necessario resolver o seguinte sistema

linear

S1(0) = 0

S ′1(0) = 0

S ′′1 (h/2) + 4S ′

1(h/2) + 13S1(h/2) = 1

S ′′1 (h) + 4S ′

1(h) + 13S1(h) = 1

substituindo os valores ja conhecidos0 0 0 1

0 0 1 0

1.130303784 3.363146642 5.567280581 11

2.938721232 5.172905726 7.134561163 11

a1

b1

c1

d1

=

0

0

1

1

.

Para determinar os parametros a2, b2, c2, d2 e necessario resolve o seguinte sistema


linear

S2(h) = S1(h)

S ′2(h) = S ′

1(h)

S ′′2 (3h/2) + 4S ′

1(3h/2) + 13S1(3h/2) = 1

S ′′2 (2h) + 4S ′

1(2h) + 13S1(2h) = 1

substituindo os valores ja conhecidos0.0231394051 0.0812022618 0.2849601057 1

0.2436067855 0.5699202114 1 0

5.616152437 7.429277250 8.701841745 11

9.353497487 10.13226121 10.26912232 11

a2

b1

c2

d2

=

0.0260060834

0.1461165184

1

1

.

E para determinar os parametros a3, b3, c3, d3 e necessario resolve o seguinte sistema

linear

S3(h) = S2(h)

S ′3(h) = S ′

2(h)

S ′′3 (3h/2) + 4S ′

3(3h/2) + 13S3(3h/2) = 1

S ′′3 (2h) + 4S ′

3(2h) + 13S3(2h) = 1

substituindo os valores ja conhecidos0.1851152409 0.324809047 0.569920211 1

0.9744271422 1.139840423 1 1

14.34165648 13.28185761 11.83640291 11

0.77152950 16.87806646 13.40368349 11

a2

b1

c2

d2

=

0.0658469906

0.1216378262

1

1

.

Resolvendo os tres sistemas lineares acima temos todos os coeficientes dos tres po-

linomios S1(t), S2(t) e S3(t)

a1 = −0, 4483601946, a2 = −0.1461795582, a3 = 0.0211777833

b1 = 0, 4480278100, b2 = 0.1444979543, b3 = −0.1627113441

c1 = 0 c2 = 0.0993745460, c3 = 0.2864665866,

d1 = 0, d2 = −0.0106627505, d3 = −0.0484863208


portanto

S1(t) = −0, 4483601946t3 + 0, 4480278100t2

S2(t) = −0.1461795582t3 + 0.1444979543t2 + 0.0993745460t− 0.0106627505

S3(t) = 0.0211777833t3 − 0.1627113441t2 + 0.2864665866t− 0.0484863208

0 0.5 1 1.5 2 2.50

0.02

0.04

0.06

0.08

0.1

0.12

solução exataspline

(a) Aproximacao e solucao exata da EDO

0 0.5 1 1.5 2 2.50

0.002

0.004

0.006

0.008

0.01

0.012

(b) Erro maximo da aproximacao

Figura 3.2: Spline cubico para aproximacao do calculo de uma EDO

A figura 3.2(a) mostra o grafico da solucao exata da EDO e o grafico da funcao S(t)

e a figura 3.2(b) mostra diferenca absoluta entre o spline e a solucao exata.

Exemplo 3.3. Resolva atraves do metodo da Sub-particao a equacao diferencial

y(5) + a4y(4) + a3y

(3) + a2y(2) + a1y

(1) + a0y = u(t)

sob as condicoes iniciais

y(0) = y(1)(0) = y(2)(0) = y(3)(0) = y(4)(0) = 0

onde a0 = 173848, 38525, a1 = 77581, 5591, a2 = 4628, 432, a3 = 595, 12, a4 = 14, 7 e u e

a funcao degrau unitario.

Suponha que a EDO e estavel. Como a EDO e linear com coeficientes constantes,

entao considere xf = ta. Agora e preciso determinar a parte real do autovalor que possui


a maior parte real da EDO, para poder determinar o tempo de acomodacao. Para isso

escreva a matriz associada a EDO em questao

A =

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

−a0 −a1 −a2 −a3 −a4

.

Sabendo que o metodo da potencia aplicado a A converge para o autovalor que possui

o maior valor absoluto de A, vamos aplica-lo a eA, para calcular a parte real do autovalor

que possui a maior parte real de A. Aplicando o metodo da potencia a eA com 20 iteracoes

descobrimos que eλ ≈ 0, 082055325534095, donde

λ ≈ −2, 500361557591399 ≈ −2, 5

e a parte real do autovalor que possui a maior parte real da EDO do problema. Sendo

assim, o Tempo de Assentamento com precisao p = 99, 5% para o nosso problema e

ta =ln(1− 0, 995)

−2, 5= − ln(0, 005)

2, 5≈ 5, 3

2, 5= 2, 12.

Vamos resolver o problema de seis formas, usando seis valores de s distintos. Tomaremos

s = 2, s = 10, s = 14, s = 18, s = 20 e s = 21 para resolver o problema. E em cada caso

precisamos resolver os sistemas

Sk(xk) = Sk−1(xk)

S(1)k (xk) = S

(1)k−1(xk)

f(Sk(xk +14h), . . . , S

(5)k (xk +

14h), xk +

14h) = r(xk +

14h)

f(Sk(xk +24h), . . . , S

(5)k (xk +

24h), xk +

24h) = r(xk +

24h)

f(Sk(xk +34h), . . . , S

(5)k (xk +

34h), xk +

34h) = r(xk +

34h)

f(Sk(xk+1), . . . , S(5)k (xk+1), xk + h) = r(xk+1)

para k = 1, 2, . . . , s.


A figura 3.3 mostra o grafico da solucao exata da EDO e do spline de aproximacao

calculado pelo metodo da Sub-particao, para os valores de s = 2, s = 10, s = 14, s = 18,

s = 20 e s = 21.

Note que o numero de pontos de maximo e mınimo locais crescem de acordo com o

valor de s, que e mais perceptıvel nas solucoes com s = 18, s = 20 e s = 21.

Como o metodo da Sub-particao descarta as condicoes iniciais da segunda derivada

adiante entao se espera que para EDOs de grau alto o spline calculado pelo metodo nao

seja uma boa aproximacao.

Pela figura 3.3 pode-se observar que o aumento do numero de divisoes do intervalo,

s, nao torna a solucao melhor, ou seja, nao adianta diminuir o passo em cada iteracao,

h, para melhorar a aproximacao atraves do metodo. Porem para o exemplo dado pode-

se observar que a solucao onde s = 2, define bem a acomodacao da solucao exata do

problema.


0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(a) solucao com s = 2

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(b) solucao com s = 10

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(c) solucao com s = 14

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(d) solucao com s = 18

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(e) solucao com s = 20

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(f) solucao com s = 21

Figura 3.3: Solucao Exata (tracejado) e spline calculado pelo metodo da Sub-particao


Outra consideracao numerica e como controlar a oscilacao da solucao limitando as

frequencias ligadas com o valor imaginario dos autovalores.

Pode-se usar uma outra adaptacao do metodo da potencia para este fim. Para evitar

que o numero de pontos de maximo e mınimo locais do spline calculado seja muito maior

que os da solucao exata, pode-se calcular o autovalor da EDO que possui a maior parte

imaginaria. Aplicando o metodo classico (algoritmo 2.1) a matriz eiA, e possıvel determi-

nar um limitante para esse valor. Com isso, podemos determinar efetivamente o numero

de intervalos necessarios para o spline. Nesse trabalho nao trataremos dessa condicao.

3.3 Metodo da Suavidade

Assim como o metodo da Sub-particao, o metodo da Suavidade sera definido de ma-

neira que possa calcular a solucao de uma EDO de ordem n de forma aproximada e tem

o objetivo que a funcao seja a mais suave possıvel, sendo nesse caso de classe Cm onde

no mınimo m = n − 1. Esse e o principal metodo desse trabalho e foi estudado mais

detalhadamente por apresentar melhores resultados que poderao ser vistos adiante. Esse

metodo e sua adaptacao serao apresentados pela primeira vez nesse trabalho.

Considere a seguinte EDO

f(y, y(1), . . . , y(n), t) = r(t)

estamos interessados em determinar uma funcao (spline) que seja uma aproximacao para

solucao exata dessa EDO, y(t).

Nesse metodo o spline devera satisfazer as condicoes iniciais e para isso seu grau deve

ser no mınimo n− 1, visto que a equacao diferencial que devera ser resolvida tem grau n.

E interessante que o spline tenha grau n− 1 para evitar calculos desnecessarios, a menos

quando se queira que a funcao seja de classe maior que Cn−1.

Da mesma forma que o metodo da Sub-particao sera preciso definir: [0, xf ], o intervalo


onde sera definido o spline, que caso a EDO seja linear com coeficientes constantes pode-se

considerar xf = ta; s, o numero total de nos do spline; h, o comprimento de cada intervalo;

xk = (k− 1)h para k = 1, 2, . . . , s as abcissas dos nos do spline; Sk para k = 1, 2, . . . , s, o

polinomio que define o spline em cada intervalo [xk, xk+1].

A diferenca entre os dois metodos esta nos sistemas de equacoes que determinam os

coeficientes dos polinomios Sk para k = 1, 2, . . . , s, enquanto o metodo da Sub-particao

exige que o spline satisfaca a EDO em diversos pontos, o metodo da Suavidade exige que

spline seja suave, ou seja, que ele possua diversas derivadas contınuas.

Para isso basta resolver os seguintes sistemas para k = 1, 2, . . . , s

Sk(xk) = Sk−1(xk)

S(1)k (xk) = S

(1)k−1(xk)

S(2)k (xk) = S

(2)k−1(xk)

......

...

S(n−1)k (xk) = S

(n−1)k−1 (xk)

f(Sk(xk+1), S(1)k (xk+1), . . . , S

(n)k , xk+1(xk+1)) = r(xk+1)

(3.4)

(Sistema de equacoes para o metodo da Suavidade)

donde S0(t) satisfaz todas as condicoes iniciais.

As n primeiras linhas do sistema 3.2 exige do spline que ele seja contınuo e que sua

derivada seja contınua no intervalo [0, xf ] e a ultima linha do sistema 3.2 exige que o spline

satisfaca a EDO no pontos xk para k = 1, 2, . . . , s. Nas n primeiras linhas do sistema 3.2

sempre temos equacoes lineares, ja na ultima linha teremos uma equacao linear se a EDO

f(y,y(1),...,y(n))(t) = r(t) for linear e nao linear caso a EDO seja nao linear.

Algoritmo do metodo da Suavidade

1. Escolha um valor para xf , caso a EDO em questao for linear com coeficientes cons-

tantes pode-se determinar que xf = ta, e calcular ta (tempo de acomodacao) usando


o Metodo da Potencia para calcular a parte real do autovalor que possui a maior

parte real da EDO.

2. Escolha o valor de s de forma conveniente e calcule o valor de h pela equacao

sh = xf .

3. Defina xk = (k − 1)h.

4. Defina o polinomio S0(x) de maneira que este satisfaca todas as condicoes iniciais.

5. Determine os coeficientes dos polinomios Sk(x), para k = 1, 2, . . . , s atraves do

sistema linear 3.4.

6. Defina o spline S(t) = Sk(t) se xk−1 < t ≤ xk.

Exemplo 3.4. Resolva usando o metodo da Suavidade a equacao diferencial linear com

coeficientes constantes

y(5) + a4y(4) + a3y

(3) + a2y(2) + a1y

(1) + a0y = u(t)


y(0) = y(1)(0) = y(2)(0) = y(3)(0) = y(4)(0) = 0

onde a0 = 173848, 38525, a1 = 77581, 5591, a2 = 4628, 432, a3 = 595, 12, a4 = 14, 7 e u e

a funcao degrau unitario. Sabendo que esta EDO e estavel.

O tempo de acomodacao dessa EDO ja foi calculado no exemplo 3.3, ta = 2, 12 = xf .

Vamos resolver esse problema de seis formas distintas cada uma com um valor de s

diferente, s = 20, s = 40, s = 60, s = 80, s = 120 e s = 160. Para cada um dos casos e

preciso resolver os sistemas


Sk(xk) = Sk−1(xk)

S(1)k (xk) = S

(1)k−1(xk)

S(2)k (xk) = S

(2)k−1(xk)

S(3)k (xk) = S

(3)k−1(xk)

S(4)k (xk) = S

(4)k−1(xk)

f(Sk(xk+1), S(1)k (xk+1), . . . , S

(5)k (xk+1), xk+1) = r(xk+1)

para k = 1, 2, . . . , s.

A figura 3.4 mostra o grafico da solucao exata da EDO (linha tracejada) e os splines

(linha cheia) calculados atraves do metodo da Suavidade, para os valores de s = 20,

s = 40, s = 60, s = 80, s = 120 e s = 160.

Note que a medida que o valor de t aumenta o spline vai piorando sua aproximacao

da solucao exata.


0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(a) solucao com s = 20

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(b) solucao com s = 40

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(c) solucao com s = 60

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(d) solucao com s = 80

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(e) solucao com s = 120

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(f) solucao com s = 160

Figura 3.4: Solucao Exata (tracejado) e spline calculado pelo metodo da Suavidade


Implementacao do metodo da Suavidade

Resolver o sistema de equacoes 3.4 para cada valor de k = 1, 2, . . . , s pode ser muito

custoso quando se quer uma melhor aproximacao tomando o valor de s alto. Porem esse

processo pede ser simplificado se for observado um padrao nesse conjunto de sistemas.

O sistema de equacoes 3.4 pode ser escrito na forma Akvk = Akvk−1

f(Sk(xk+1), S′k(xk+1), . . . , S

(n)k (xk+1), xk+1) = r(xk+1)

(3.5)

onde vk representa a solucao do sistema de equacoes 3.4 (coeficientes do spline) e Ak e

uma matriz de ordem n× (n+ 1) para k = 1, 2, . . . , s.

Definindo Cjn como

Cjn =

n

j

=n!

j!(n− j)!,

pode-se verificar que Ak em 3.5 tem a forma

Ak =

xnk xn−1

k xn−2k . . . x2

k xk 1

nxn−1k (n− 1)xn−2

k (n− 2)xn−3k . . . 2xk 1 0

2C2nx

n−2k 2C2

n−1xn−3k 2C2

n−2xn−4k . . . 2 0 0

......

.... . .

......

...

(n− 2)!C(n−2)n x3

k (n− 2)!C(n−2)n−1 x2

k (n− 2)!C(n−2)n−2 xk . . . 0 0 0

(n− 1)!C(n−1)n x2

k (n− 1)!C(n−1)n−1 xk (n− 2)! . . . 0 0 0

n!xk (n− 1)! 0 . . . 0 0 0

.

O seguinte resultado sera usado para simplificar o algoritmo.

Lema 3.1. As matrizes M de ordem n × (n + 1) e N(h) de ordem (n + 1) × (n + 1)

definidas por


M =

0 0 0 . . . 0 0 1

0 0 0 . . . 0 1 0

0 0 0 . . . 2 0 0...

......

. . ....

......

0 0 (n− 2)! . . . 0 0 0

0 (n− 1)! 0 . . . 0 0 0

(3.6)

e

N(h) =

1 0 0 . . . 0 0 0

nh 1 0 . . . 0 0 0

C2nh

2 (n− 1)h 1 . . . 0 0 0

C3nh

3 C2n−1h

2 (n− 2)h . . . 0 0 0...

......

. . ....

......

Cn−2n hn−2 Cn−3

n−1hn−3 Cn−4

n−2hn−4 . . . 1 0 0

Cn−1n hn−1 Cn−2

n−1hn−2 Cn−3

n−2hn−3 . . . 2h 1 0

hn hn−1 hn−2 . . . h2 h 1

(3.7)

satisfazem

Ak = MN(h)−(k−1)

onde a matriz Ak esta definida em 3.5.

Demonstracao. Sejam M (em 3.6) e N(h) (em 3.7). A matriz N(h) e inversıvel (sua

inversa satisfaz N(h)−1 = N(−h)).

Rapidamente pode-se verificar que A1 = M , visto que x1 = 0.

Pode-se verificar ainda que Ak+1 = AkN(h)−1 e assim concluir que

Ak = MN(h)−(k−1). (3.8)

Assim, a partir do lema 3.1, defina M e N = N(h)−1. Como no sistema 3.5 temos

Akvk = Akvk−1


logo o lema 3.1 nos da

MNk−1(vk − vk−1) = 0

assim

Nk−1(vk − vk−1) = αke1

onde

e1 =[1 0 . . . 0 0

]te αk ∈ R

visto que αke1 ∈ Ker(M).

Dessa forma

vk = vk−1 + αkN−(k−1)e1

onde αk ainda nao foi calculado.

Agora falta apenas determinar o valor de αk considerando a condicao de interpolacao

(segunda linha do sistema 3.5) que e possivelmente nao linear.

Assim para aplicar o metodo da Suavidade e necessario determinar a variavel αk

resolvendo a equacao

f(Sk(xk+1), S′k(xk+1), . . . , S

(n)k (xk+1), xk+1) = r(xk+1)

que para simplifica-la sera denotada por

f(vk, xk+1) = r(xk+1)

ou

f(vk−1 + αkuk−1, xk+1) = r(xk+1). (3.9)

onde uk−1 representa a primeira coluna de N−(k−1), ou seja, uk−1 = N−(k−1)e1.

Quando a EDO for linear a equacao 3.9 tambem e linear e podemos simplificar a

solucao do sistema 3.4 a produtos entre matrizes e vetores que veremos a seguir. Quando

a EDO e nao linear a equacao 3.9 tambem e nao linear e portanto sera mais difıcil resolve-

la tendo que recorrer a metodos numericos de resolucao de equacoes nao lineares com uma

variavel.


Alem disso, quando a EDO for nao linear a equacao podera ter mais de uma solucao

e nao se sabe se a solucao e a ideal ou nao.

Sendo assim trataremos aqui apenas o caso linear, ou seja,

y(n) + pn−1(t)yn−1 + . . .+ p1(t)y

′ + p0(t)y = r(t)

sendo pi(t) constantes ou nao, o processo para determinar uma aproximacao para y(t)

pelo metodo da Suavidade pode ser simplificado.

Definindo Sk(t) (polinomio de grau n) como

Sk(t) = vkntn + vk(n−1)t

(n−1) + . . .+ vk1t+ vk0

onde vkj ∈ R para todo k = 1, 2, . . . , s e todo j = 0, 1, . . . , n, para aplicar o metodo da

Suavidade, e preciso resolver a ultima equacao do sistema 3.4 dada por

S(n)k + pn−1(t)S

n−1k + . . .+ p1(t)S

′k + p0(t)Sk = r(t)

para k = 1, 2, . . . , s.

Agora definindo os vetores

fxk+1=

p0(xk+1)xnk+1 + p1(xk+1)nx

n−1k+1 + . . .+ pn−1(xk+1)n!xk+1 + n!

p0(xk+1)xn−1k+1 + p1(xk+1)(n− 1)xn−2

k+1 + . . .+ pn−1(xk+1)(n− 1)!...

p0(xk+1)xk+1 + p1(xk+1)

p0(xk+1)

t

e

vk =[vkn vk(n−1) . . . vk1 vk0

]tde modo que

fk+1vk = r(xk+1)

temos

fk+1(vk−1 + αkuk−1) = r(xk+1)


logo

fk+1vk−1 + αkfk+1uk−1 = r(xk+1)

e portanto

αk =r(xk+1)− fk+1vk−1

fk+1uk−1

.

Com isso podemos descrever o algoritmo completo.


Algoritmo 3.1: metodo da Suavidade Implementado Para EDO’s Lineares

Entrada: p0(t), p1(t), . . . , pn−1(t), r(t), y0, y(1)0 , . . . , y

(n−1)0 , xf ,s.

Saıda: S(t)

inıcio

h← xf/s;

para j=0 ate n-1 faca

v1j ← y(j)0 ;

fim para

v1n ← (r(0)− y0 − y′0 − . . .− y(n)0 )/n!;

para i=0 ate n faca

para j=0 ate n faca

se j > i entao Nij ← 0;

senao Nij ← Ci−jn−j+1(−h)i−j;

fim para

fim para

u1 ← [Nj0, Nj1, . . . , Njn]t;

para k=2 ate s faca

xk+1 = kh;

define fk+1;

αk = (r(xk+1)− fk+1vk−1)/(fk+1uk−1) ; /* "condic~ao crıtica" */

vk = vk−1 + αkuk−1;

uk = Nuk−1;

fim para

fim


O unico lugar no algoritmo 3.1 onde entram os coeficientes da EDO linear (que podem

ser ou nao constantes) e na linha marcada como ”condicao crıtica”.

Outra maneira de simplificar a resolucao do sistema 3.4 para EDO lineares e definindo

o vetor

gk+1 =[n! pn−1(xk+1) . . . p1(xk+1) p0(xk+1)

]te como Ak

fk+1

vk =

Akvk−1

r(xk+1)

e Ak

fk+1

=

M

gk+1

Nk−1

assim a solucao do sistema 3.4 e dado por

vk = N−(k−1)

M

gk+1

−1 Akvk−1

r(xk+1)

.

Essa ultima forma pode parecer mais simples a primeira vista, mas o custo computaci-

onal e um pouco maior visto que e necessario realizar produtos do tipo matriz-matriz en-

quanto a forma anterior realiza no maximo operacoes do tipo matriz-vetor a cada iteracao

(Golub e Van Loan [5]).

Uma das qualidades dos algoritmos spline, conforme os exemplos ja discutidos e que a

solucao aproximada converge no regime permanente ao valor assintotico correto da solucao

exata.

Exemplo 3.5. Conside a EDO linear com coefiecientes constantes

y′′ + 5y′ + 4y = u(t) (3.10)

com condicoes iniciais y(0) = y′(0) = 0.

Definindo as variaveis de entrada

p0(t) = 4, p1(t) = 5, r(t) = u(t), y(0) = y′(0) = 0, xf = 10 e s = 100


de acordo com o algoritmo 3.1 temos

h =xf

s= 0.02, v1 =

[1 0 0

]te u1 =

[1 −nh C2

nh2

]t.

Os resultados das iteracoes para alguns valores de k de 1 ate 100 estao presentes na

tabela 3.3.

k xk ak vk

1 0 - (0.5, 0, 0)

2 0.2 -0.54 (−0.4 10−1, 0.216,−0.216 10−1)

3 0.4 -0.368 10−1 (−0.768 10−1, 0.24544,−0.27488 10−1)...

......

...

49 9.6 0.225640 10−5 (−0.1047488 10−4, 0.22475882 10−3, 0.24878336724)

50 9.8 0.185649 10−5 (−0.861839 10−5, 0.18837153 10−3, 0.24896166497)

51 10.0 0.152746 10−5 (−0.709093 10−5, 0.15782228 10−3, 0.24911441123)...

......

...

98 19.4 0.15 10−9 (−0.73 10−9, 0.3034 10−7, 0.24999968774)

99 19.6 0.13 10−9 (−0.60 10−9, 0.2521 10−7, 0.24999973806)

100 19.8 0.10 10−9 (−0.50 10−9, 0.2094 10−7, 0.24999978032)

Tabela 3.1: Resolucao do exemplo 3.5 pelo metodo da Suavidade com s = 100 e xf = 20

Pode-se observar, nesse exemplo, que valor absoluto de ak diminui a medida que o

valor de k aumenta levando os ultimos polinomios do spline a serem bem proximos.

No capıtulo 1 vimos que, para a EDO da equacao 3.10

lim y(t) =1

4= 0.25

assim nao e surpreendente que em uma aproximacao para y definida a partir de polinomios

do tipo

Sk = vk2t2 + vk1t+ vk0

as constantes vk2 e vk1 convirjam para 0 e vk0 convirja para 1/4.


A solucao da EDO da equacao 3.10 e a aproximacao (spline) calculado pelo metodo

da Suavidade Implementado, assim como a diferenca absoluta dessa duas solucoes pode

ser visto na figura 3.5.

0 5 10 15 200

0.05

0.1

0.15

0.2

0.25

t

Am

plitu

de

(a) Solucao e Spline (s = 100)

0 5 10 15 20−5

0

5

10

15

20x 10

−3

t

Err

o

(b) Erro Absoluto (s = 100)

Figura 3.5: Solucao Exata (tracejado) e spline calculado pelo metodo da Suavidade Im-

plementado

3.4 metodo da Suavidade Adaptado

Considerando que o spline calculado a partir do metodo da Suavidade aproxima-se

bem da solucao exata de uma EDO para os valores de t iniciais, mas que a medida

que o valor de t aumenta essa aproximacao torna-se cada vez pior, enquanto o spline

calculado a partir do metodo da Sub-particao, para valores de s pequenos aproxima-se

bem melhor que o spline calculado a partir do metodo da Suavidade para os valores de t

mais proximos do tempo de acomodacao, podemos considerar um metodo que seja uma

adaptacao do metodo da Suavidade onde a medida que o valor de t aumenta seja exigida

menos suavidade do spline. Alem disso, em uma EDO estavel quando t → ∞ teremos

y(n), y(n−1), . . . , y(2), y(1) → 0.

Para determinar um spline que a medida que o valor de t aumenta seu grau diminua,

primeiro precisa-se determinar como sera feita essa diminuicao e depois em como dividir


o intervalo [0, xf ], onde para EDO lineares com coeficientes constantes pode-se tomar

xf = ta. Para a diminuicao do grau do spline propomos que a cada novo intervalo o

spline diminua em 1 o seu grau ate chegar ao mınimo de grau 2, visto que e interessante

que o spline seja contınuo e sua derivada seja contınua. Portanto se a EDO tiver grau n,

precisamos dividir o intervalo [0, xf ] em n− 1 intervalos. Sendo assim dada uma equacao

diferencial ordinaria linear de grau n

f(y, y(1), . . . , y(n), t) = r(t)

e definido o valor final xf , onde o spline sera definido, dividiremos o intervalo [0, xf ] em

n − 1 intervalos de mesmo comprimento, que chamaremos de blocos e denotaremos por

bk para k = 1, 2, . . . , n − 1, donde cada bloco tera comprimento h dado pela equacao

(n−1)h = ta. E cada bloco bk dividiremos em sk intervalos, para k = 1, 2, . . . , n−1, para

que em cada bloco tenhamos sk polinomios de grau n− k+1. Sendo assim a medida que

o valor de k aumenta o grau do polinomio sera cada vez menor e chegara ao valor mınimo

de 2 quando k = n− 1.

Em seguida e necessario que se defina os valores de sk, para k = 1, 2, . . . , n−1, e como

a ideia inicial e que a exigencia sobre a suavidade do spline diminua a medida que o valor

de t aumenta, o numero de spline em cada intervalo deve satisfazer

s1 > s2 > . . . > sn−1.

Apos a definicao do numero de nos em cada spline basta aplicar o metodo da Suavidade

a cada um dos blocos bk para k = 1, 2, . . . , n − 1, onde em cada bloco bk sera calculado

um spline de grau n− k + 1.

Por exemplo, para resolver uma EDO linear com coeficientes constantes estavel de

grau 6 pelo metodo da Suavidade Adaptado, primeiro iremos determinar o tempo de

acomodacao, ta. Depois dividir o intervalo [0, ta] em 5 intervalos de comprimento h = ta/5,

sao eles: b1, b2, b3, b4 e b5. Com isso podemos definir os valores para s1, s2, s3, s4 e s5,

considerando a ordem decrescente podemos escolher os valores s1 = 16, s2 = 8, s3 =

4, s4 = 2 e s5 = 1. O spline do bloco b1 tera grau 6, o spline do bloco b2 tera grau 5, o


spline do bloco b3 tera grau 4, o spline do bloco b4 tera grau 3 e o spline do bloco b5 tera

grau 2.

Figura 3.6: Exemplo de como aplicar o metodo da Suavidade Adaptado a uma EDO de

grau 6

Como no bloco bk teremos sk intervalos, entao cada um desses intervalos tera compri-

mento hk = h/sk, para k = 1, 2, 3, 4, 5. A figura 3.6 mostra como fica o intervalo [0, ta]

apos a sua divisao em blocos e a divisao dos blocos em intervalos.

Algoritmo do metodo da Suavidade Adaptado

1. Defina intervalo [0, xf ], onde sera definido o spline. Para EDOs lineares com coe-

ficientes constantes pode-se usar o Metodo da Potencia para calcular a parte real

do autovalor que possui a maior parte real da EDO e com isso calcular o tempo de

acomodacao, e definir xf = ta.

2. Divida o intervalo [0, xf ] em n−1 intervalos (blocos) e determine o valor de h, donde

h = xf/(n− 1).

3. Determine o numero de intervalos, sk, em cada um dos blocos bk, de forma que

s1 > s2 > . . . > sn−1.

4. Determine o comprimento de cada um dos intervalos dos blocos, hk = h/sk.

5. Determine o spline de grau n−k+1 em cada um dos blocos bk para k = 1, 2, . . . , n−1

pelo metodo da Suavidade.


6. Defina o spline S(t) como sendo cada um dos splines calculados nos blocos bk para

k = 1, 2, . . . , n− 1.

Exemplo 3.6. Resolva usando o metodo da Suavidade Adaptado a equacao diferencial

y(5) + a4y(4) + a3y

(3) + a2y(2) + a1y

(1) + a0y = u(t)


y(0) = y(1)(0) = y(2)(0) = y(3)(0) = y(4)(0) = 0

onde a0 = 173848, 38525, a1 = 77581, 5591, a2 = 4628, 432, a3 = 595, 12, a4 = 14, 7 e u e

a funcao degrau unitario. Sabendo que esta EDO e estavel.

O tempo de acomodacao dessa EDO ja foi calculado no exemplo 3.3, ta = 2, 12.

Vamos resolver esse problema de seis formas distintas cada uma com um conjunto de

valores {s1, s2, s3, s4} diferentes. Para cada um dos casos e preciso resolver os sistemas

Si,k(xk) = Si,k−1(xk)

S(1)i,k (xk) = S

(1)i,k−1(xk)

S(2)i,k (xk) = S

(2)i,k−1(xk)

S(3)i,k (xk) = S

(3)i,k−1(xk)

S(4)i,k (xk) = S

(4)i,k−1(xk)

f(Si,k(xk+1), S(1)i,k (xk+1), . . . , S

(5)i,k (xk+1), xk+1) = r(xk+1)

para k = 1, 2, . . . , si, donde i = 1, 2, 3, 4.

A figura 3.7 mostra o grafico da solucao exata da EDO (linha tracejada) e os splines

(linha cheia) calculados atraves do metodo da Suavidade, para os quatro conjuntos de

valores distintos, {s1, s2, s3, s4}.


0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(a) s1 = 8, s2 = 4, s3 = 2 e s4 = 1

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(b) s1 = 16, s2 = 8, s3 = 4 e s4 = 2

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(c) s1 = 32, s2 = 16, s3 = 8 e s4 = 4

0 0.5 1 1.5 20

1

2

3

4

5

6

x 10−6

(d) s1 = 64, s2 = 32, s3 = 16 e s4 = 8

Figura 3.7: Solucao Exata (tracejado) e spline calculado pelometodo da Suavidade Adap-

tado


3.5 Resultados Numericos

Nessa secao serao feitas comparacao entre os metodos: metodo da Suavidade, metodo

de Euler e calculo dos autovalores para EDO’s lineares sendo que este ultimo sera consi-

derado apenas caso a EDO seja linear com coeficientes constantes. Nas comparacoes nao

serao considerados os outros dois metodos (da Sub-particao e da Suavidade Adaptado),

pois estes nao convergem para a solucao exata a medida que aumentamos o tamanho da

particao, visto que nao consideram todas as condicoes de suavidade.

A forma da solucao exata de EDO’s lineares com coeficientes constantes e conhecida

e depende exclusivamente dos autovalores da EDO e das condicoes iniciais. Ao resolver

uma EDO por um desses metodos numericos pode-se comparar com a solucao dada pelo

calculo dos autovalores em tempo e aproximacao da solucao exata. Nessa secao serao

feitas essas comparacoes.

Para o calculo dos autovalores usamos o Matlab que utiliza dois metodos: decom-

posicao de Cholesky (ver, por exemplo, Watkins [13]), para matrizes simetricas, e decom-

posicao QZ para autovalores generalizados no caso geral. A decomposicao QZ tambem e

conhecida como decomposicao generalizada de Schur.

Um estudo mais detalhado sobre metodo de Euler para resolucao de uma EDO pode

ser encontrado em Boyce e DiPrima [2]. O algoritmo 3.2 descreve os passos necessarios

para resolver uma EDO numericamente por esse metodo.

Sera denotado At(n) para representar o tempo medio para resolver uma EDO linear

com coeficientes constantes atraves do calculo dos seus autovalores, e Et(n, s) e St(n, s)

para representar, respectivamente, o tempo medio para resolver uma EDO linear com

coeficientes constantes de grau n com s iteracoes pelo metodo de Euler e pelo metodo da

Suavidade.


Algoritmo 3.2: Metodo de Euler

Entrada: f(y, y′, . . . , y(n−1), t), x0, xf , y0, y(1)0 , . . . , y

(n−1)0 , s.

Saıda: y0, y1, . . . , ys+1

inıcio

dx← (xf − x0)/s;

para i=0 ate s faca

para j=0 ate n-1 faca

y(j)i+1 ← y

(j)i + fj(yi, y

′i, . . . , y

(n−1)i , xk)dx;

fim para

xi+1 ← xk−1 + dx;

fim para

fim

O tempo relativo para resolver uma EDO linear com coeficientes constantes atraves

do metodo de Euler sera definido como

Etr(n, s) =Et(n, s)

At(n)

e atraves do metodo da Suavidade sera definido como

Str(n, s) =St(n, s)

At(n).

Se S(t) e a solucao aproximada de uma EDO pelo metodo da Suavidade, dados h e s,

definimos como erro relativo do metodo da Suavidade como

Ser =

∫ sh

0|S(t)− y(t)| dt∫ sh

0|y(t)| dt

onde y(t) e a solucao exata da EDO.

Se yk e a solucao aproximada da EDO no ponto xk para k = 1, 2, . . . , s pelo metodo

de Euler, definimos como erro relativo do metodo de Euler como

Eer =

∑sk=0 |yk+1 − yk|h/2 dt∫ sh

0|y(t)| dt

onde y(t) e a solucao exata da EDO.


Exemplo 3.7. Considere a seguinte EDO de grau 4,

y(4) + 8y(3) + 25y(2) + 46y(1) + 40y = u(t)

que possui autovalores

λ1 = −4, λ2 = −2, λ3 = −1 + 2i, λ4 = −1− 2i

A EDO sera resolvida no intervalo [0, 6] para diversos valores de s pelo metodo de

Euler e pelo metodo da Suavidade. A tabela 3.5 mostra o tempo computacional e o erro

relativo das solucoes pelo metodo de Euler e pelo metodo da Suavidade em relacao ao

tempo computacional para determinar a solucao exata calculando os autovalores.

s Etr Str Eer Ser Et/St Eer/Ser

24 0.8281 2.3754 0.0209 0.0037 0.3486 5.5846

38 1.0660 2.7813 0.0089 0.0024 0.3833 3.6483

50 1.2898 3.1731 0.0060 0.0019 0.4065 3.1830

Tabela 3.2: Comparacao entre as solucoes do exemplo 3.7 pelo metodo de Euler e pelo

metodo da Suavidade com o calculo dos autovalores

De acordo com a tabela 3.5 pode-se concluir que calculo aproximado pelo metodo de

Euler e pelo metodo da Suavidade nao sao competitivos com o calculo exato da solucao

da EDO do exemplo 3.7.

Para comparar o metodo de Euler com o metodo da Suavidade ao resolver a EDO do

exemplo 3.7 e interessante considerar solucoes que se tenha o tempo computacional ou o

erro relativo proximos.

Nessa situacao se ao resolver a EDO do exemplo 3.7 pelo metodo da Suavidade for

escolhido um valor de s menor que 120 e ao resolver essa EDO pelo metodo de Euler for

escolhido um valor de s apropriado, entao a solucao pelo metodo de Euler tera erro e

tempo computacional menor.


Suavidade Euler Et/St Eer/Ser

s = 30 s = 85 0.7375 1.0094

s = 30 s = 122 0.9910 0.6608

s = 120 s = 287 0.9944 1.0002

s = 1000 s = 1950 1.0081 1.1571

s = 1000 s = 2254 1.1628 1.0001


metodo da Suavidade

Analogamente, se ao resolver a EDO do exemplo 3.7 pelo metodo de Euler escolher um

valor de s maior que 287 e ao resolver essa EDO pelo metodo da Suavidade for escolhido

um valor de s apropriado, entao a solucao pelo metodo da Suavidade tera erro e tempo

computacional menor.

Os valores de 120 e 287 sao aproximados e foram escolhidos de acordo com os valores

de Et e St (tempos medios) que foram calculados resolvendo o problema diversas vezes e

calculando a media do tempo.

0 1 2 3 4 5 60

0.005

0.01

0.015

0.02

0.025

0.03

Am

plitu

de

t

EulerSuavidadeExata

(a) Solucoes Aproximadas e Exata

0 1 2 3 4 5 60

1

2

3

4

5

6

7

8x 10

−4

Err

o A

bsol

uto

t

EulerSuavidade

(b) Erro Absoluto

Figura 3.8: Solucao da EDO do exemplo 3.7 pelo metodo de Euler com s = 175 e da

Suavidade com s = 70

A figura 3.8 mostra o grafico das solucoes do metodo de Euler e do metodo da Sua-

vidade para valores de s onde o tempo computacional e o erro relativo entre as solucoes


sao praticamente iguais.

Exemplo 3.8. Considere a seguinte EDO linear

y′′ + ty′ + 4y = 2 t cos(2t)

com condicoes iniciais y0 = 0 e y′0 = 2 que possui solucao exata igual a

y(t) = sen(2t).

A EDO sera resolvida no intervalo [0, 6] para diversos valores de s pelo metodo de Euler

e pelo metodo da Suavidade. A tabela 3.5 mostra as relacoes entre o tempo computacional

e o erro relativo das solucoes aproximadas calculadas pelo metodo de Euler e pelo metodo

da Suavidade.

s Et/St Eer Ser Eer/Ser

20 0.3287 3.0852 0.5853 5.2715

50 0.3873 0.9339 0.2450 3.8116

90 0.4338 0.4861 0.1379 3.5251


metodo da Suavidade e pelo calculo dos autovalores

As solucoes desse problema pelo metodo da Suavidade e pelo metodo de Euler com

valores de s = 20 e s = 90 podem ser observado na figura 3.9 assim como o erro absoluto

em cada caso. As figuras (a) e (c) apresentam o grafico da solucao exata da EDO e os

graficos com as solucoes aproximadas calculadas pelo metodo de Euler e pelo metodo da

Suavidade.

Na figura 3.9 pode-se observar que nos graficos apresentados o erro absoluto das

solucoes calculadas pelo metodo da Suavidade e menor que o erro das solucoes calcu-

ladas pelo metodo de Euler, porem, e importante ressaltar que o tempo computacional

gasto para calcular as solucoes pelo Metodo da Euler foi, em cada um dos casos, menor


0 1 2 3 4 5 6−2.5

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2A

mpl

itude

t

EulerSuavidadeExata

(a) Solucao com s = 20

0 1 2 3 4 5 60

0.2

0.4

0.6

0.8

1

1.2

1.4

Err

o A

bsol

uto

t

EulerSuavidade

(b) Erro Absoluto com s = 20

0 1 2 3 4 5 6−1.5

−1

−0.5

0

0.5

1

1.5

Am

plitu

de

t

EulerSuavidadeExata

(c) Solucao como s = 90

0 1 2 3 4 5 60

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Err

o A

bsol

uto

t

EulerSuavidade

(d) Erro Absoluto com s = 90

Figura 3.9: Solucao da EDO do exemplo 3.8 pelo metodo de Euler e da Suavidade

que o tempo gasto para calcular as solucoes pelo metodo da Suavidade e isso pode ser

observado na tabela 3.5.

Para comparar as solucoes entre os dois Metodos (Euler e Suavidade) e interessante

que o tempo computacional ou o erro relativo sejam iguais, ou bem proximos, e para isso

basta tomar valores de s convenientes em cada caso.

Resolvendo o problema dado no exemplo 3.8 para diversos valores de s pelo metodo

da Suavidade e pelo metodo de Euler foi observado que ao resolver o problema pelo

metodo da Suavidade com o valor de s maior de 52 podem ser escolhidos valores de s

para o metodo de Euler que mostram que a solucao pelo metodo da Suavidade possui


Suavidade Euler Et/St Eer/Ser

s = 20 s = 76 0.6980 0.9976

s = 20 s = 120 0.9933 0.6110

s = 30 s = 108 0.8148 0.9968

s = 30 s = 139 0.9915 0.7637

s = 52 s = 178 1.0000 1.0018

s = 300 s = 628 1.0021 1.5659

s = 300 s = 979 1.5440 1.0009

s = 1000 s = 1910 1.0019 1.7001

s = 1000 s = 3242 1.6967 1.0003


metodo da Suavidade

erro e tempo computacional menor, analisando a razao entre os tempos computacionais

e a razao entre os erros relativos de acordo com a tabela 3.5.

Para exemplificar, considere resolver o problema do exemplo 3.8 pelo metodo da Sua-

vidade com s = 300 e para comparar resolver o problema, tambem, pelo metodo de Euler

com s = 628 e depois com s = 979. No primeiro caso, os tempos computacionais para a

resolucao sao bem proximos, mas o erro relativo da aproximacao calculado pelo metodo

de Euler e bem maior. No segundo caso, os erros relativos sao bem proximos, porem o

tempo computacional da resolucao pelo metodo de Euler e bem maior. Nos dois casos o

metodo da Suavidade foi mais eficiente.

Analogamente ao escolher um valor de s menor que 52 ao resolver o problema pelo

metodo da Suavidade pode-se observar que as solucoes pelo metodo de Euler serao me-

lhores em erro e tempo computacional.

Conclusao

Para determinar uma solucao aproximada para a EDO o metodo da Sub-particao nao

se mostrou um bom metodo, visto que este desconsidera as condicoes iniciais para EDO’s

de grau maiores que 2 e diminuir a particao, aumentando o valor de s, nao melhora a

aproximacao. Porem considerando EDO’s lineares com coeficientes constantes estaveis

este metodo mostrou-se bom para estudar a acomodacao da solucao exata.

Para resolver EDOs lineares com coeficientes constantes de grau baixo (considerando

grau menor que 20) o calculo numerico dos autovalores mostrou-se superior a qualquer

metodo desenvolvido nesse trabalho considerando tempo computacional e erro relativo

para equacoes de grau baixo; porem para equacoes de grau alto pode-se considerar re-

solver a EDO atraves do metodo da Suavidade, ou pelo metodo de Euler, visto que o

tempo computacional para calcular os autovalores de uma matriz cresce muito mais rapi-

damente que o tempo para realizar as operacoes elementares de uma iteracao do metodo

da Suavidade e do metodo de Euler.

Alem disso, e uma boa escolha optar por resolver uma EDO linear com coeficientes

constantes pelo calculo numerico dos autovalores quando e necessario estudar todo o

regime transitorio da EDO, ademais quando esta possui um tempo de acomodacao muito

alto devido a margem pequena de estabilidade. Em contrapartida quando o tempo de

acomodacao for muito baixo ou quando for necessario estudar a solucao apenas em um

intervalo pequeno os metodos da Suavidade e de Euler serao boas opcoes visto que o

tempo computacional podera ser bem menor.

82


De toda forma o metodo da Suavidade mostra-se bastante util visto que resolve bem

nao so EDO lineares com coeficientes constantes como tambem EDO lineares com coefici-

entes variaveis e o metodo da Suavidade resolve essas EDO’s com praticamente a mesma

dificuldade de uma equacao com coeficientes constantes, a menos de avaliacoes das funcoes

que definem a EDO em cada iteracao.

Ao comparar o metodo da Suavidade, apresentado nesse trabalho, com o metodo de

Euler, um metodo classico para resolucao de EDO’s numericamente foi visto que depen-

dendo da necessidade de uma aproximacao melhor, tomando o valor de s alto, o metodo

da Suavidade sera mais eficiente que o metodo de Euler, porem quando se deseja uma

aproximacao simples, quando o valor de s e baixo, o metodo de Euler e uma opcao melhor.

Nessa dissertacao foi feito um trabalho preliminar com o metodo da Suavidade, reali-

zando uma descricao do metodo, uma implementacao dele e analises numericas compara-

tivas com outros metodos. Porem alguns resultados importantes nao foram determinados

como: calculo de numero de operacoes e avaliacoes em funcao do grau da EDO e do numero

de passos; um limitante superior do erro relativo acumulado, como funcao da malha; e

uma comparacao rigorosa com os outros metodos (incluindo metodos de Runge-Kutta de

maior ordem), que leve em consideracao simultaneamente o erro e a complexidade (tempo

computacional). Esses resultados tornariam o metodo mais significativo e sao possıveis

caminhos para a continuidade do estudo do tema.

Apendice A

Propriedades Basicas da

Transformada de Laplace

A Transformada de Laplace de uma funcao satisfaz as seguintes propriedades:

p1 : L[a1y1 + a2y2] = a1L[y1] + a2L[y2], (linearidade)

p2 : L[y(k)] = skL[y]− sk−1y(0)− sk−2y(1)(0)− ...− sy(k−2)(0)− y(k−1)(0)

p3 : L[e−at] =

1

s+ a

onde a, a1, a2 ∈ C e y, y1 e y2 sao funcoes definidas para todo t ≥ 0 e y(i) denota a

i-esima derivada de y.

A linearidade da Tranformada de Laplace e herdada da linearidade da integral de uma

funcao.

A propriedade p2 pode ser demonstrada com o princıpio da inducao finita, basta aplicar

a integracao por partes. Para k = 1, temos

L[y(1)] =

∫ ∞

0

e−sty(1)(t)dt = e−sty(t)∣∣∞0−∫ ∞

0

−se−sty(t)dt

= −y(0) + s

∫ ∞

0

e−sty(t)dt = sL[y]− y(0).

84

Apendice A. Propriedades Basicas da Transformada de Laplace 85

Agora suponha que a propriedade e valida para k = n− 1, ou seja,

L[y(n−1)] = sn−1L[y]− sn−2y(0)− sn−3y(1)(0)− ...− sy(n−3)(0)− y(n−2)(0)

assim

L[y(n)] =

∫ ∞

0

e−sty(n)(t)dt = e−sty(n−1)(t)∣∣∞0−∫ ∞

0

−se−sty(n−1)(t)dt

= −y(n−1)(0) + s

∫ ∞

0

e−sty(n−1)(t)dt = sL[y(n−1)]− y(n−1)(0)

= snL[y]− sn−1y(0)− sn−2y(1)(0)− ...− sy(n−2)(0)− y(n−1)(0).

A propriedade p3 pode ser verificada calculando L(e−at) diretamente da definicao de

Transformada de Laplace, ou seja, se Re(s+ a) > 0

L[e−at] =

∫ ∞

0

e−ste−atdt =

∫ ∞

0

e−(s+a)tdt =1

s+ a.

Para mais detalhes sobre a solucao de EDOs Lineares com coeficientes constantes ver

Madureira [8].

Apendice B

Estabilidade Assintotica de EDOL’s

com Coeficientes Constantes

Proposicao B.1. Se a parte real de todos os autovalores da EDO 1.7 e negativa, entao

quando t→∞ sua solucao converge para 1/a0.

Demonstracao. Se a parte real de todos os autovalores da EDO 1.7 e negativa entao

cada parcela Akeλkt e uma funcao que converge para zero, pois escrevendo λk = ck + dki

onde ck, dk ∈ R, temos, pela identidade de Euler, que

Akeλkt = Ake

ckt(cos(dkt) + isen(dkt))

mas Ak(cos(dkt) + isen(dkt)) e uma funcao limitada e eckt converge para zero, entao

limt→∞

Akeλkt = 0

e como Akeλkt converge para zero, para k = 1, 2, . . . , n entao, pela proposicao 1.4, temos

limt→∞

y(t) = A0.

Proposicao B.2. Se a parte real de um dos autovalores da EDO 1.7 e positiva entao

quando t→∞ sua solucao diverge.

86

Apendice B. Estabilidade Assintotica de EDOL’s com Coeficientes Constantes 87

Demonstracao. Seja λ1 = c1 + d1i, λ2 = c2 + d2i, . . ., λn = cn + dni os autovalores da

EDO, organizados de modo que c1 ≥ c2 ≥ . . . ≥ cn e se c1 > 0 (hipotese), de acordo com

o corolario 1.1 podemos escrever a solucao dessa EDO na forma

y(t) = A0 +n∑

k=1

eckt(ukcos(dkt) + vksen(dkt)).

suponha que λ1 = λ2, ou seja, que λ1 /∈ R (a demonstracao para λ1 = λ2 e similar). Pela

proposicao 1.2 temos u1 = u2 e v1 = −v2 logo

y(t) = A0 + ec1t(v1cos(d1t) + v1sen(d1t) +n∑

k=3

e(ck−c1)t(ukcos(dkt) + vksen(dkt)))

e como e(ck−c1)t → 0 quando t → ∞ e (ukcos(dkt) + vksen(dkt)) e limitado para todo

k = 3, 4, . . . , n, entao quando t→∞

e(ck−c1)t(ukcos(dkt) + vksen(dkt))→ 0, para k = 3, 4, . . . , n

logo

limt→∞

y(t) = A0 + limt→∞

ec1t(v1cos(d1t) + v1sen(d1t))

e portanto y(t) diverge.

Para mais detalhes sobre a solucao de EDOs Lineares com coeficientes constantes ver

Boyce e DiPrima [2].

Apendice C

Programacao no MATLAB

Todos os resultados em que houve necessidade de calculo computacional nesse trabalho

(incluindo todos os graficos apresentados) foram determinados a partir de programacao

usando o software Matlab v7.8.0.347 (R2009a) 64 bit em um sistema windows 7 (seven).

Um referencial pratico para desenvolver a programacao, os graficos e realizar alguns

calculos desse trabalho pode ser encontrado em Ogata [11]. Em Mathews [9] pode-se en-

contrar resultados sobre a parte teorica de alguns desses calculos como metodos numericos

para a resolucao de equacoes diferenciais e o metodo da potencia.

Segue a programacao de alguns algoritmos apresentados nesse trabalho na linguagem

MATLAB:

1. Resolucao do problema de resposta a degrau;

2. Metodo da Potencia (Classico);

3. Metodo da Potencia (Adaptado);

4. Metodo de Euler;

5. Metodo da Suavidade.

88

Apendice C. Programacao no MATLAB 89

1 - Resolucao do problema de resposta a degrau

% ---------------------------------------------------------------

% Codigo para resolver uma EDO linear com coeficientes constantes

a = [1 6 538];

u = [0 1];

TF = tf(u,a);

[y,t] = step(TF);

plot(t,y,’b’);

xlabel(’t’);

ylabel(’amplitude’);

% ---------------------------------------------------------------


2 - Metodo da Potencia (Classico)

% ---------------------------------------------------------------

function av = metPot(A,nInt)

% Metodo da Potencia aplicada a Matriz A

if (nargin < 2) nInt = 50; end

n = length(A);

x = rand(n,1)+i*rand(n,1);

y = x;

for k=1:1:nInt

xaux = x(:,k);

yaux = y(:,k);

yprox = (1/max(xaux))*xaux;

xprox = A*yprox;

x = [x xprox];

y = [y yprox];

end

av = real((xprox’*yprox)/(yprox’*yprox));

end

% ---------------------------------------------------------------


3 - Metodo da Potencia (Adaptado)

% ---------------------------------------------------------------

function lk = metPotAdap(A,nInt)

% Metodo da Potencia Adaptado aplicada a Matriz A

if (nargin < 2) nInt = 50; end

n = length(A);

u = zeros(nInt,n);

N = zeros(nInt,1);

u(1,:) = rand(n,1);

u(1,:) = [0.1808 0.2456 0.7964 0.8558];

N(1) = 1/norm(u(1,:));

u(1,:) = u(1,:)*N(1);

for i=2:1:nInt

u(i,:) = A*u(i-1,:)’;

N(i) = 1/norm(u(i,:));

u(i,:) = u(i,:)*N(i);

end

z0 = u(nInt,:)’;

z1 = A*z0;

z2 = A*z1;

R = [z0(1) z1(1); z0(2) z1(2)];

solve = R\[z2(1);z2(2)];

a = solve(1);

b = solve(2);

parteReal = b/2;

parteImaginaria = sqrt(b^2 + 4*a)/2;

lk = parteReal + parteImaginaria;

end

% ---------------------------------------------------------------


4 - Metodo de Euler

% ---------------------------------------------------------------

function [xvalues, yvalues] = euler(f,x0,xn,y0,s)

% Metodo de Euler para resolver EDO’s

dx = (xn-x0)/s;

x = x0:dx:xn-dx;

n = length(y0);

y = zeros(s,n);

y(1,:) = y0;

for k=1:1:s-1

y(k+1,:)= y(k,:) + f(y(k,:),x(k))*dx;

end

xvalues = x’;

yvalues = y(:,1);

end

% ---------------------------------------------------------------


5 - Metodo da Suavidade

% ---------------------------------------------------------------

function V = metodoSuavidade(f,r,x0,xf,y0,s)

% O Metodo da Suavidade retorna um vetor V com os coeficientes do Spline

h = (xf-x0)/s;

x = linspace(x0,xf,s+1);

n = length(y0);

v(1,:) = [(r(0) - sum(y0))/factorial(n) y0(n:-1:1)];

N = zeros(n+1,n+1);

for i=1:1:n+1

for j=1:1:i

N(i,j) = (-h)^(i-j)*nchoosek(n-j+1,i-j);

end

end

u(1,:) = N(:,1)’;

for k = 2:1:s

F = f(x(k));

v(k,:) = v(k-1,:) + ((r(x(k)) - F*v(k-1,:)’)/(F*u(k-1,:)’))*u(k-1,:);

u(k,:) = N*u(k-1,:)’;

end

V = v;

end

% ---------------------------------------------------------------

Referencias Bibliograficas

[1] AHLBERG, J. H., NILSON, E. N. e WALSH, J. L. The Theory of Splines and

Their Applications, vol. 38. Editora Academic Press Inc, New York, 1967.

[2] BOYCE, William E. e DiPRIMA, Richard C. Elementary differential equations

and boundary value problems, 7 ed. Ed. Wiley, 2000.

[3] BURDEN, Richard L. e FAIRES, J. Douglas. Analise Numerica, 1 ed. Editora

Thomson, 2003

[4] COHEN, Nir e LEWKOWICZ, Izchak, The Lyapunov Order for Real Matrices.

Linear Algebra Appl. 2009.

[5] GOLUB, Gene H. e VAN LOAN, Charles F. Matrix Computations, 3 ed. Editora

The Johns Hopkins University Press, Baltimore, 1996.

[6] HORN, Roger A. e JOHNSON, Charles R. Matrix Analysis, Cambridge University

Press, 1985.

[7] LATHI, Bhagawandas P. Sinais e Sistemas Lineares, 2 ed. Editora Artmed, Sao

Paulo, 2004.

[8] MADUREIRA, Luısa. Problemas de Equacoes Diferenciais Ordinarias e

Transformadas de Laplace, vol. 3. Editora FEUP Edicoes, Porto, 2010.

[9] MATHEWS, John H. Numerical Methods: MATLAB Programs, 1995.

94

Referencias Bibliograficas 95

[10] OGATA, Katsuhiko. Engenharia de Controle Moderno., 3 ed. Rio de Janeiro,

LTC, 2000.

[11] OGATA, Katsuhiko.Matlab For Control Enginners, 1. ed. Editora Prentice Hall,

2008.

[12] POOLE, David. Linear Algebra: A Modern Introduction, 3 ed. Editora Brooks

Cole, 2010.

[13] WATKINS, David S. Fundamentals of Matrix Computations, John Wiley &

Sons, NY, 1991.

Documents

Métodos Numéricos para Resolu¸cão de Equa¸cões ... · Resumo Neste trabalho desenlvolvemos um método de resolu¸cão de problemas de valor inicial ... Um artigo recente,