Métodos_de_Integração_Numérica

Prévia do material em texto

Teoria de Controle e suas aplicações 
 
197 
15 Métodos de integração numérica 
A integração numérica ganhou importância quando grandes sistemas não 
lineares passaram a descrever sistemas físicos de grande porte. Entre eles podemos citar 
as grandes redes elétricas modeladas por suas equações algebrico-diferenciais, com 
fenômenos não lineares representados. Nestes casos, a solução tem que 
obrigatoriamente ser obtida através de simulação no domínio do tempo, pois só 
podemos empregar a solução no domínio da frequência em sistemas lineares, onde é 
aplicada a transformada de Laplace e são empregadas regras da Álgebra Linear para 
fazer a análise do comportamento dinâmico do sistema, sem a necessidade de 
simulações no domínio do tempo. 
15.1 Método de Euler simples 
O método de Euler simples, também conhecido por Euler explícito, é o mais 
simples entre todos os métodos de integração numérica e também o que mais produz 
erros de 2ª discretização. É chamado de explícito, pelo fato de no cálculo do valor da 
variável de saída, em um determinado instante, depender apenas dos valores da entrada 
e saída no instante imediatamente anterior. 
Inicialmente vamos observar como devemos proceder para resolver uma 
equação diferencial de 1ª ordem qualquer, dada por: 
) ( ty,f
dt
dy
 
A solução analítica dessa equação é encontrada aplicando-se a integração. Isto é: 
Teoria de Controle e suas aplicações 
 
198 
Cdtty,fty   ) ()( 
Por esse motivo que os métodos numéricos utilizados para resolver esse tipo de 
equação são chamados de métodos de integração. 
Podemos observar que na sua solução, a constante C depende de condições de 
contorno, para que esta possa ser determinada. Caso contrário, a solução será dada por 
uma família de curvas, cuja solução é calculada pela integral indefinida da função 
) ( ty,f . 
Para a solução numérica dessa equação, sempre será necessário fornecer as 
condições iniciais, obtidas através das condições de contorno. As condições iniciais 
servem para começar o processo de integração, cujo cálculo é feito pela área 
compreendida entre a curva e o eixo das abscissas, onde geralmente é o eixo que 
representa o tempo. 
Nos métodos numéricos de integração é feita uma aproximação dos termos 
diferenciais para diferenças, discretizando o problema. Os erros relacionados com essa 
aproximação são chamados de erros de 2ª discretização, uma vez que os de 1ª 
discretização estão relacionados com os erros que o computador comete na 
representação dos números reais e complexos. 
Então, com as aproximações mencionadas, a integral utilizada na solução da 
equação diferencial se transforma em um somatório de parcelas discretizadas no tempo, 
cuja parcela de discretização passa de diferencial de tempo (dt) para uma diferença de 
tempo (t). Essa diferença de tempo também é muito conhecida como passo de 
integração (h). 
Teoria de Controle e suas aplicações 
 
199 
Com isso a solução de ) ( ty,f
dt
dy
 passa de 
maxt
t
dtty,fty
1
) ()( para 



maxt
tt
tty,fty
1
) ()( , que é dada pela soma acumulada das áreas de retângulos de altura 
) ( ty,f e largura t, conforme está ilustrado na Figura 15.1. 
 
Figura 15.1 – Integração pelo Método de Euler Simples 
 
Como o problema foi discretizado e amostrado no tempo, a solução numérica do 
problema é dada por uma sequência de valores sob forma de tabela onde podemos 
atribuir um índice p para cada um deles, indicando o respectivo passo de integração, 
como segue: 
 
 
Teoria de Controle e suas aplicações 
 
200 
Tabela 15.1 – Evolução da solução do método de Euler simples 
p tempo y 
1 t1 y1 
2 t2 = t1 + h y2 = y1 + h f(y1, t1) 
   
p tp+1 = tp + h yp+1 = yp + h f(yp, tp) 
   
np hnthtt p
n
max
p )1(11   yn+1 = yn + h f(yn, tn) 
 
Facilmente podemos observar que o número total de passos de integração (np) é 
dado por: 
1
1



h
tt
n maxp 
É muito comum se fazer a observação gráfica do resultado, uma vez que é muito 
difícil analisar o comportamento da solução através de valores numéricos. Para isso, 
basta ligarmos os pontos obtidos por seguimentos de reta. 
Uma forma gráfica para vizualizar como é feito o cálculo da área abaixo da 
curva em um passo de integração pelo método de Euler simples é ilustrada na Figura 
15.2. Podemos verificar que a referida área é obtida pelo cálculo da área do retângulo e 
que o erro cometido por esse método é bem grosseiro, cujo valor está representado por 
. No gráfico da esquerda, que representa a solução, as respectivas áreas do gráfico da 
direita se transformam em ordenadas, pois há um processo de integração. 
Teoria de Controle e suas aplicações 
 
201 
 
Figura 15.2 – Método de Euler Simples 
 
Podemos apresentar o que foi mostrado graficamente na Figura 15.2 sob forma 
de algoritmo a ser implementado em computador digital, da seguinte maneira: 













 Façado Fim
 
) ( 
 2 1Faça 
 :iniciais Valores
1
1
11
htt
t,yfhyy
h
t
,,,p
t,y
pp
pppp
max
 
Para generalizar o problema, vamos agora considerar um sistema de equações 
diferenciais de 1ª ordem descrito por: 













) (
) (
) (
21
212
2
211
1
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
nn
n
n
n




 
Teoria de Controle e suas aplicações 
 
202 
Podemos estender o algoritmo da solução apresentada para uma única equação 
diferencial pelo método de Euler simples para a seguinte forma: 























 Façado Fim
 
 Façado Fim
) ( 
 , 2 1Faça 
 
 2 1Faça 
 :iniciais Valores
1
21
1
111
2
1
1
htt
t,y,,y,yhfyy
n,,k
h
t
,,,p
t,y,,y,y
pp
pp
n
pp
k
p
k
p
k
max
n




 
Vamos considerar agora uma outra forma para a equação diferencial de 1ª 
ordem, dada por: 
u=a y+y  
Sendo y a variável de estado e u a entrada do sistema. Integrando a equação no 
intervalo de t a (t+t), isto é, em um passo de integração, obtemos a seguinte expressão: 

 tt
t
tt
t
tt
t
u dt = a y dt + dy 
A sua solução é dada por: 
ttu=ttya+tytty  )( )( )()( 
E, portanto, chegamos a: 
 )( )() 1()(
21
tuttyta=tty
CC
  
Mais uma vez, podemos observar que o método é explícito porque o valor de 
)( tty  só depende dos valores de y(t) e y(t), que já são conhecidos. 
Podemos fazer uma representação da solução do sistema em notação 
discretizada, amostrada em intervalos h, da seguinte forma: 
Teoria de Controle e suas aplicações 
 
203 
ppp uCy=Cy 21
1  e 1 h=tt pp  
15.2 Método de Euler implícito 
O método de Euler implícto é muito utilizado nos casos em que possíveis 
oscilações numéricas podem atrapalhar o cálculo do processo de integração. O grande 
problema da sua utilização está no fato de além de amortecer as oscilações numéricas 
indesejáveis, ele também amortece as oscilações reais, mascarando uma possível 
instabilidade do sistema. Portanto, a sua utilização requer uma análise mais aprofundada 
do sistema a ser resolvido. 
Este método é muito parecido com o método de Euler simples. A diferença 
consiste basicamente no cálculo da área entre a curva que representa a função 
) ( ty,f
dt
dy
 e o eixo dos tempos. Essa área é calculada pelo retângulo com base h e 
altura ) ( ht,yf *  . Como a altura desse retângulo depende do valor da variável y, no 
instante t + h, isto é, no mesmo instante de cálculo da variável, há a necessidade de se 
aplicar um processo iterativo até que a diferença entre o valor de y calculado e o 
respectivo valor anterior seja menor que uma dada tolerância. Por esse motivo é que 
esse método de integração faz parte da família dos métodos implícitos. 
Uma forma gráfica para vizualizar como é feito o cálculo da área abaixo da 
curva em um passo de integração pelo método de Euler implícito é apresentada na 
Figura 15.3.Podemos verificar que a referida área é obtida pelo cálculo da área do 
retângulo, cujo valor do erro está representado por . No gráfico da esquerda, que 
representa a solução, as respectivas áreas do gráfico da direita se transformam em 
ordenadas, pois há um processo de integração. 
Teoria de Controle e suas aplicações 
 
204 
 
Figura 15.3 – Método de Euler Implícito 
Então, podemos esboçar um algoritmo correspondente à utilização do método de 
Euler implícito, a seguir: 


















































 Façado Fim
 
enquanto Faça do Fim
 
1faça então E0 Se 
) ( 
0 
1 enquantoFaça 
 
1 
) ( 
 2 1Faça 
 :iniciais Valores
1
1
1
1
1
11
htt
yy
conv,tolerância
y
yy
conv
t,yhfyy
conv
conv
conv
t,yhfyy
h
t
,,,p
,ty
pp
pant
p
antp
pantpp
pppant
max
 
Observamos que, no início de cada passo de integração, é feito o cálculo de uma 
estimativa para o valor de y no final do intervalo, para que o processo de convergência 
Teoria de Controle e suas aplicações 
 
205 
possa ser iniciado. Nesse algoritmo fizemos esse cálculo tomando como base a derivada 
da variável y no início do intervalo. Porém, nem sempre essa é a melhor estratégia para 
se obter o menor número médio de iterações por passo, quando se cobre todo o tempo 
de simulação. 
Sistema de equações diferenciais de 1ª ordem: 













) (
) (
) (
21
212
2
211
1
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
nn
n
n
n




 
Algoritmo da solução pelo método de Euler implícito: 















































































 Façado Fim
 
enquanto Faça do Fim
 Façado Fim
 
1faça então E0 Se 
) ( 
 , 2 1Faça 
 
0 
1 enquantoFaça 
 
1 
 Façado Fim
) ( 
 , 2 1Faça 
 
 2 1Faça 
 :iniciais Valores
1
1
1
1
21
1
21
111
2
1
1
htt
yy
conv,tolerância
y
yy
conv
t,y,,y,yhfyy
n,,k
conv
conv
conv
t,y,,y,yhfyy
n,,k
h
t
,,,p
t,y,,y,y
pp
p
k
ant
k
p
k
ant
k
p
k
pant
n
antant
k
p
k
p
k
pp
n
pp
k
p
k
ant
k
max
n






 
Teoria de Controle e suas aplicações 
 
206 
Vamos considerar agora uma outra forma para a equação diferencial de 1ª 
ordem: 
u=a y+y  
Sendo y a variável de estado e u a entrada do sistema. Integrando a equação no 
intervalo de t a (t + t), isto é, em um passo de integração, obtemos a seguinte 
expressão: 

 tt
t
tt
t
tt
t
u dt = a y dt + dy 
A sua solução é dada por: 
tttu=tttya+tytty  )( )( )()( 
E, portanto, podemos chegar a: 
tttut=yttyta  )( )()() 1( 
Ou, então: 
)(
 1
)(
 1
1
)(
21
ttu
ta
t
ty
ta
=tty
CC

















 
Podemos observar mais uma vez que o método é implícito, porque o valor de 
)( tty  não depende só do valor de y(t), que é conhecido, mas também depende do 
valor de )( ttu  , que não é conhecido. Esta característica é uma das razões da 
necessidade da solução ser iterativa. Para diminuir o número de iterações deste método 
de integração implícito, podemos aplicar a extrapolação quadrática de u(t), que é mais 
eficiente em relação às demais. Provaremos mais adiante que a sua expressão é dada por 
)2()]()([3)( ttuttutu=ttu  , desde que o passo de integração seja constante 
nos instantes )( e ),( ),2( ttttttt  . 
Teoria de Controle e suas aplicações 
 
207 
Podemos fazer uma representação da solução do sistema em notação 
discretizada, amostrada em intervalos h, da seguinte forma: 
1
21
1   ppp uCy=Cy e 1 h=tt pp  
A Figura 15.4 mostra o polinômio de segundo grau a ser utilizado na 
extrapolação quadrática para o cálculo do valor inicial de )( ttu  , em função dos 
valores conhecidos das entradas nos três instantes anteriores a )( tt  , isto é, 
)2( e )( ),( ttuttutu  . 
 
Figura 15.4 – Extrapolação Quadrática 
 
A expressão geral da função quadrática no tempo é dada por: 
CBtAttu  2)( 
Então, podemos utilizar os três pontos nos instantes anteriores a )( tt  para 
formar o sistema de equações lineares e obtermos os valores dos coeficientes A, B e C. 








 )2()2()2(
)()()(
 )(
2
3
2
2
2
1
CttBttAttuu
CttBttAttuu
CBtAttuu
 
Na forma matricial, temos: 
Teoria de Controle e suas aplicações 
 
208 



































 
 
 
 
1)2()2(
1)()(
1
3
2
1
2
2
2
u
u
u
C
B
A
tttt
tttt
tt
 
 
Do sistema de equações anterior podemos obter as expressões de A, B e C, 
utilizando a regra de Cramer, como segue: 
2
321
2
2
2
3
2
1
)(2
2
 
1)2()2(
1)()(
1
 
 
1)2(
1)(
1
 
t
uuu
tttt
tttt
tt
ttu
ttu
tu
A







 
2
321
2
2
2
3
2
2
2
1
2
)(2
 )2( )(4 )32(
 
1)2()2(
1)()(
1
 
 
1)2(
1)(
1
 
t
uttuttutt
tttt
tttt
tt
utt
utt
ut
B







 
2
3
2
2
2
1
22
2
2
2
3
2
2
2
1
2
)(2
 )( )4(2 ])(23[
 
1)2()2(
1)()(
1
 
 
)2()2(
)()( 
t
utttutttutttt
tttt
tttt
tt
utttt
utttt
utt
C







 
Com os valores de A, B e C determinados podemos obter a expressão de 
)( ttu  , substituindo-os na equação geral, como segue: 
Teoria de Controle e suas aplicações 
 
209 
2
3
2
2
2
1
22
2
321
2
2
321
)(2
 )( )4(2 ])(23[
 
)(
)(2
 )2( )(4 )32(
 
)(
)(2
2
)(
t
utttutttutttt
tt
t
uttuttutt
tt
t
uuu
ttu









 
Então, chegamos na equação final: 
)2()]()([3)(3)( 321 ttuttutuuuu=ttu  
Podemos observar o baixo esforço computacional para se fazer a extrapolação 
quadrática, isto é, utiliza-se apenas uma operação de multiplicação e duas adições. O 
ganho no desempenho computacional compensa fazer essas operações no início de cada 
passo de integração, pois diminui bastante o número global de iterações. Na aplicação 
da extrapolação quadrática é exigido que se tenha 3 valores consecutivos já calculados, 
com isso no início ou em algum instante que haja descontinuidade devido a possíveis 
perturbações aplicadas ao sistema, há a necessidade de utilizar alguma extratégia na 
obtenção desses 3 valores. Uma delas pode ser não se fazer nenhum tipo de 
extrapolação, isto é, pode-se partir o processo de integração utilizando o mesmo valor 
do instante anterior. Isso não onera o tempo global de simulação, pois a extrapolação 
não será usada somente no instante inicial e nos instantes em que aconteçam 
descontinuidades no sistema. 
15.3 Método de Euler modificado 
No método de Euler modificado a área de cada passo de integração é calculada 
pela área do trapézio com base h e alturas ) ( ty,f e ) ( ht,yf *  , nos instantes t e (t+h), 
respectivamente, conforme mostrado na Figura 15.5. Também é um método implícito, 
Teoria de Controle e suas aplicações 
 
210 
pois a altura ao final do intervalo depende do valor de y, que não é conhecido. Então, no 
processo iterativo de otenção de y, o seu valor inicial é calculado pelo método de Euler 
simples. A convergência do processo iterativo é obtida no momento que a diferença 
entre o valor de y calculado e o respectivo valor anterior seja menor que uma dada 
tolerância. 
 
Figura 15.5 – Método de Euler Modificado 
Podemos esboçar um algoritmo correspondente à utilização do método de Euler 
modificado, como segue: 
Teoria de Controle e suas aplicações 
 
211 





























































 Façado Fim
 
enquanto Faça do Fim
2
 
1faça então E0 Se 
)(
2
 
)( 
 
0 
1 enquantoFaça 
 
1 
) ( 
 2 1Faça 
 :iniciais Valores
1
1
1
1
11
htt
BA
A
conv,tolerância
y
yy
conv
BA
h
yy
hhA,tyfB
hAyy
conv
conv
conv
t,yfA
h
t
,,,p
,ty
pp
p
antp
pp
pp
pant
pp
max
 
Observamos que, no início de cada passo de integração, é feito o cálculo de uma 
estimativa para o valor de y no final do intervalo, para que o processo de convergência 
possa ser iniciado. Nesse algoritmo, para representar o método de Euler modificado, 
esse cálculo é feito tomando como base a derivada da variável y no início do intervalo. 
Porém, nem sempre essa é a melhor estratégia para se obter o menor número médio de 
iterações por passo, quando se cobre todo o tempo de simulação. 
Para o problema mais geral, vamos considerar agora um sistema de equações 
diferenciais de 1ª ordem descrito por: 
Teoria de Controle e suas aplicações 
 
212 













) (
) (
) (
21
212
2
211
1
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
t,y,,y,yf
dt
dy
nn
n
n
n




 
Podemos estender o algoritmo da solução apresentada para uma única equação 
diferencial pelo método de Euler modificado para a seguinte forma: 





































































































 Façado Fim
 
enquanto Faça do Fim
 Façado Fim
2
 
1faça então E0 Se 
)(
2
 
) ( 
 
 , 2 1Faça 
 
0 
1 enquantoFaça 
 
1 
 Façado Fim
) ( 
 , 2 1Faça 
 
 2 1Faça 
 :iniciais Valores
1
1
1
1
21
21
111
2
1
1
htt
BA
A
conv,tolerância
y
yy
conv
BA
h
yy
t,y,,y,yfB
hAyy
n,,k
conv
conv
conv
t,y,,y,yfA
n,,k
h
t
,,,p
t,y,,y,y
pp
kk
k
p
k
ant
k
p
k
kk
p
k
p
k
pant
n
antant
kk
k
p
k
ant
k
pp
n
pp
kk
max
n






 
Teoria de Controle e suas aplicações 
 
213 
Vamos considerar agora uma outra forma para a equação diferencial de 1ª 
ordem, a seguir: 
u=a y+y  
Sendo y a variável de estado e u a entrada do sistema. Integrando a equação no 
intervalo de t a (t+t), isto é, em um passo de integração, obtemos a seguinte expressão: 

 tt
t
tt
t
tt
t
u dt = a y dt + dy 
A sua solução é dada por: 
)]()([
2
 )]()([ 
2
 )()( tuttu
t
=tytty
t
a+tytty 



 
E, portanto chegamos a: 
)(
2
1
2 )(
2
1
2 )(
2
1
2
1
 )(
221
ttu
t
a
t
+tu
t
a 
t
+t y
t
a 
t
a 
=tty
CCC











































 
O termo b1 é dependente apenas dos valores históricos y(t) e u(t), enquanto que o 
termo b2 é uma função de u(t+t). 
Podemos, então, fazer uma representação da solução do sistema em notação 
discretizada, amostrada em intervalos h, da seguinte forma: 
1
221
1   pppp uCuCy=Cy e 1 h=tt pp  
Observamos que o método é implícito, porque o valor de )( tty  não depende 
só dos valores de y(t) e u(t), que são conhecidos, mas também depende do valor de 
)( ttu  , que não é conhecido. Esta característica é uma das razões da necessidade da 
solução ser iterativa. Para diminuir o número de iterações deste método de integração 
implícito, podemos aplicar a extrapolação quadrática de u(t), que é mais eficiente em 
Teoria de Controle e suas aplicações 
 
214 
relação às demais. Quando isso é feito, o método de Euler modificado passa a ser 
chamado de método trapezoidal implícito. Se não for feito o processo iterativo e nem 
algum tipo de extrapolação, esse método de 2ª ordem passa a ser chamado de método de 
Runge-Kutta de 2ª ordem. A ordem dos métodos de integração está associada com o 
erro global produzido pela soma dos erros de cada passo de integração. Uma análise 
superficial pode ser feita considerando que uma função possa ser desenvolvida em série 
de Taylor, como segue: 





!3
)(
)(
2
)(
)()()()(
32 t
tg
t
tgttgtgttg '''''' 
Se associarmos 
dt
dy
 a )()( tftg'  , então )(1 ttgy p  , )( pp tgy  , 
)()(
2
2
p'p'' tftg
dt
yd
 , )()(
3
3
p''p''' tftg
dt
yd
 e assim por diante. Substituindo na 
equação desenvolvida em série de Taylor, passamos a ter: 
 )(
!3
)(
!2
)(
32
1 p''p'ppp tf
h
tf
h
thfyy 
Com essa equação, podemos ver que os métodos de 1ª ordem só utilizam os 
dois primeiros termos da série. Com isso, o erro cometido em cada passo de integração 
é dado pela soma dos termos seguintes. Podemos observar que os termos posteriores 
sempre serão menores que os anteriores, pois o coeficiente de cada termo é formado por 
uma função potência de h no numerador e uma função fatorial do índice do termo no 
denominador. Isso faz com que o coeficiente de cada termo tenda para zero a medida 
que o índice do termo cresce. Então, podemos dizer que o erro cometido em um passo 
de integração é da ordem de h2, pois está associado com o maior termo da série do seu 
cálculo. Como o erro global é dado pela soma dos erros de cada passo de integração, 
Teoria de Controle e suas aplicações 
 
215 
dizemos que ele é da ordem de h para os métodos de 1ª ordem. Para os métodos de 2ª 
ordem, o erro global é da ordem de h2, e assim por diante. 
15.4 Método de Runge-Kutta de 4ª ordem 
Para o cálculo do valor da variável de estado no passo de integração seguinte, o 
método de Runge-Kutta de 4ª Ordem utiliza a média ponderada dos valores das 
derivadas da variável de estado em quatro pontos, como apresentado na Figura 15.6. 
 
Figura 15.6 – Método de Runge-Kutta de 4ª ordem 
Teoria de Controle e suas aplicações 
 
216 
O primeiro é ponto ) ,( pp ty , isto é, no início do intervalo. O segundo é o ponto 
)2/ ,2/( hthAy pp  , localizado no meio do intervalo. O terceiro, também localizado 
no meio do intervalo, é o ponto )2/ ,2/( hthBy pp  . O quarto é ponto 
) ,( hthCy pp  , está no final do intervalo. Inicialmente vamos observar como 
devemos proceder para resolver uma equação diferencial de 1ª ordem qualquer, dada 
por: 
)( y,tf
dt
dy
 


























htt
DCBA
h
yy
ht,hCyfD
h
t,B
h
yfC
h
A,t
h
yfB
,tyfA
h
t
,,,p
,ty
pp
pp
pp
pp
pp
pp
max
1
1
11
)22(
6
)(
)
22
(
)
22
(
)(
21
 :iniciais Valores

 
Seja agora um sistema de equações diferenciais de 1ª ordem descrito por: 













)(
)(
)(
21
212
2
211
1
,tx,,,xxf
dt
dx
,tx,,,xxf
dt
dx
,tx,,,xxf
dt
dx
nn
n
n
n




 
 
A solução pelo método de Runge-Kutta de 4ª ordem é dada a seguir: 
Teoria de Controle e suas aplicações 
 
217 
 





































































htt
DCBA
h
xx
DCBA
h
xx
DCBA
h
xx
ht,hCx,,hCx,hCxfD
ht,hCx,,hCx,hCxfD
ht,hCx,,hCx,hCxfD
h
t,B
h
x,,B
h
x,B
h
xfC
h
t,B
h
x,,B
h
x,B
h
xfC
h
t,B
h
x,,B
h
x,B
h
xfC
h
t,A
h
x,,A
h
x,A
h
xfB
h
t,A
h
x,,A
h
x,A
h
xfB
h
t,A
h
x,,A
h
x,A
h
xfB
t,x,,x,xfA
t,x,,x,xfA
t,x,,x,xfA
h
t
,,,p
t,x,,x,x
pp
nnnn
p
n
p
n
pp
pp
p
n
p
n
pp
nn
p
n
p
n
pp
p
n
p
n
pp
p
n
p
n
pp
nn
p
n
p
n
pp
p
n
p
n
pp
p
n
p
n
pp
nn
p
n
p
n
pp
p
n
p
n
pp
pp
n
pp
nn
pp
n
pp
pp
n
pp
max
n
1
1
22222
1
2
11111
1
1
2211
221122
221111
2211
221122
221111
2211
221122
221111
21
2122
2111
111
2
1
1
)22(
6
)22(
6
)22(
6
) (
) (
) (
)
2
 
2
 
2
 
2
(
)
2
 
2
 
2
 
2
(
)
2
 
2
 
2
 
2
(
)
2
 
2
 
2
 
2
(
)
2
 
2
 
2
 
2
(
)
2
 
2
 
2
 
2
(
) (
) (
) (
 2 1:iniciais Valores



















 
 
 
Teoria de Controle e suas aplicações 
 
218 
15.5 Estimativa para o valor do passo de integração 
O passo de integração (h), como já foi dito anteriormente, é a parcela de 
discretização no tempo que é feito para se resolver numericamente um sistema 
envolvendo equações diferenciais. O seu valor deve ser escolhido adequadamente de 
modo a não causar oscilações numéricas no processo de cálculo, quando este é maior do 
que um certo valor. Valores muito pequenos, além de proporcionar um esforço 
computacional muito grande, também podem ocasionar problemas numéricos, pois o 
computador comete erros de 1ª discretização na representação dos número reais, 
portanto, em termos práticos, o valor de h não deve ser menor do que 100 vezes o valor 
do  da máquina. O valor ideal para h pode ser obtido, de acordo com o método de 
integração aplicado, de tal forma que não provoque erro excessivo no cálculo da área 
abaixo da curva dada por )( y,tf
dt
dy
 , acarretando que o segmento de reta obtido da 
ligação de 2 pontos consecutivos da solução numérica possa praticamente ser 
aproximado à curva entre esses pontos. 
O problema consiste em ter uma maneira prática de se medir as aproximações 
mencionadas. É importante mencionar que esta só pode ser utilizada em sistemas 
lineares onde os seus autovalores possam ser calculados. Em sistemas não lineares não 
há alternativa para isso, isto é, a única forma de se obter o melhor valor para h é com 
um processo de tentativas, começando com um determinado valor e avaliando a 
resposta no tempo, em seguida se reduz o seu valor e novamente avalia-se a resposta no 
tempo, comparando-se as duas. Caso essas respostas sejam diferentes, prossegue-se 
diminuindo o valor de h até que as respostas sejam praticamente coincidentes, então, o 
seu valor ideal será o mediatamente superior, pois atenderá aos requisitos mencionados. 
Teoria de Controle e suas aplicações 
 
219 
Para uma melhor visualização do problema, vamos considerar que a resposta de 
um sistema linear tenha um par de autovalores complexos conjugados 1 =  + j e 2 
=  – j. Associado a esses autovalores está a resposta no domínio do tempo dada por 
tjtj eKeKty )(2
)(
1)(
   ou, em outra forma mais amigável por 
) ()(   tsenKety t . A transformação de uma forma para a outra pode ser obtida 
através do método das frações parciais aplicado na solução de equações diferenciais 
homogêneas. Pela segunda forma podemos traçar o gráfico apresentado na Figura 15.7, 
onde podem ser observados a constante de tempo 

1
T  e o período de oscilação 

 2 . Esses parâmetros são fundamentais para que possamos estimar o passo de 
integração, cujo valor deve ser um submúltiplo do menor deles, devido às razões já 
apresentadas. Portanto, podemos obter a expressão 
n
,Tmin
h
)( 
 , com o valor de n 
sendo dependente do método de integração a ser aplicado. Valores obtidos com a 
prática indicam n = 100 para o método de Euler simples, n = 20 para o método de Euler 
Modificado com apenas uma correção e n = 10 para o método de Runge-Kutta de 4ª 
ordem. A redução do valor de n e, consequentemente o aumento do valor de h, pode ser 
feita em função da precisão do cálculo da área empregado pelo método de integração. 
Quanto maior for o valor do passo de integração, melhor será o desempemho 
computacional, pois teremos um menor número de passos a serem calculados para 
cobrir o mesmo tempo máximo de simulação. 
A função min(T,) fornece o menor valor entre todas as constantes de tempo e 
todos os períodos de oscilação do sistema em estudo (que podem ter vários estados 
Teoria de Controle e suas aplicações 
 
220 
representados) deve ser utilizada para garantir que os segmentos de reta formados pelos 
pontos consecutivos da solução numérica formem a curva de resposta com uma boa 
aproximação. 
De uma maneira geral a resposta no domínio do tempo de um sistema linear é 
dada por tn
tt neKeKeKty   21 21)( , se não tiver multiplicidade de autovalores, 
caso contrário fica na forma 
t
n
t
m
tm
m
mneKeKetKtKKty  
   21 1
1
21 )()( , considerando m sendo o 
fator de multiplicidade do autovalor 1. 
 
Figura 15.7 – Constante de Tempo e Período de Oscilação 
 
O cálculo da constante de tempo pode ser feito através da função tan(), que é 
obtida pela divisão de K pelo valor de T. Também a função tan() pode ser calculada 
pela derivada da envoltória da função y(t) no instante inicial. 
Teoria de Controle e suas aplicações 
 
221 
 
  KeK
T
K
t
t 
0
 )tan( 
Com isso, podemos concluir que 

1
T . 
Por definição a constante de tempo é o valor do tempo em que a tangente à curva 
da envoltória da função y(t) no instante inicial leva para atingir o valor final de regime 
permanente y(∞). Tomando como base essa definição, podemos concluir que o valor de 
 tem que ser negativo, caso contrário não existe regime permanente, isto é, o valor 
final y(∞) tende para infinito. Podemos dizer também, pelas mesmas razões, que para 
um sistema ser estável, os seus autovalores devem estar posicionados no semiplano da 
esquerda do plano complexo. 
A definição de período de oscilação é dada pelo tempo que um valor máximo da 
função y(t) leva para atingir o valor máximo seguinte, completando um ciclo. Como 
sabemos que f 2  e que 
f
1
 , com  sendo a pulsação em rad/s, f a frequência 
em Hz e  o período de oscilação em segundos, podemos obter a relação 

 2 . 
Associado a um par de autovalores complexos conjugados está uma equação 
diferencial homogênea de 2ª ordem, que podemos escrevê-la da seguinte forma: 
02
2
2
  y
dt
dy
dt
yd
nn  
Cujo polinômio característico é 022  nn  e suas raízes são: 


 jj n
nnn 



 

 1
2
442 2
222
 
Teoria de Controle e suas aplicações 
 
222 
 
onde: 
n  
n
21 






 







22222222222 1
cos
nnn
n
nn
n
)(
 
 
A Figura 15.8 mostra em uma forma gráfica a definição do coeficiente de 
amortecimento  , relacionado com o autovalor  =  + j . 
 
 
Figura 15.8 – Coeficiente de amortecimento 
 
Então, por definição: 
ntoamortecime de ecoeficient 
0)para a (igual oscilação de natural frequência  n 
 
 
Teoria de Controle e suas aplicações 
 
223 
15.6 Estimativa para o valor do tempo máximo de simulação 
Aproveitando o que foi dito no item anterior e observando a Figura 15.7, 
podemos ver que o tempo máximo de simulação (tmax) só depende da parcela 
relacionada com o amortecimento da resposta no domínio do tempo. Nesse tipo de 
problema só estamos interessados em observar o regime transitório da resposta, portanto 
só devemos fazer os cálculos até o instante que o valor da variável de estado esteja 
praticamente constante, estabelecendo um determinado valor para uma banda morta ., 
na qual seus valores permanecem dentro. Pelo exposto, podemos estabelecer uma 
estimativa de tmax como um múltiplo da maior contante de tempo do sistema. O valor 5 é 
muito usado para esse múltiplo, que na envoltória da curva fornece um valor de 
aproximadamente 0,00674 do valor inicial, que corresponde a fazermos t = 5T, sendo 
calculado por e-5T/T = e-5. 
Então, podemos dizer que )max( 5 Ttmax  . 
A função max(T) fornece o maior valor entre todas as constantes de tempo do sistema 
em estudo, que podem ter vários estados representados. Esta função deve ser utilizada 
para que realmente todo o sistema tenha atingido o regime permanente, pois a sua 
resposta global é dada pela soma de cada resposta individual, por estarmos 
considerando um sistema linear e portanto, é válida a aplicação do teorema da 
superposição. 
15.7 Transformação de uma equação diferencial em um sistema 
Para que possamos resolver numericamente uma equação diferencialhomogênea 
de um grau n qualquer temos que tranformá-la em um sistema de equações diferenciais 
Teoria de Controle e suas aplicações 
 
224 
de 1ª ordem. Uma das formas mais usadas para se obter essa transformação é criando 
(n-1) variáveis auxiliares, sendo cada uma delas dada pela derivada em cascata da 
antecessora. 
uxa
dt
dx
a
dt
xd
a
dt
xd
a
n
n
nn
n
n  

 10
1
11
1
1
1
1  
n
n
n
n
n
nn
n
dt
xd
dt
dx
dt
xd
dt
xd
dt
dx
x
dt
xd
dt
dx
x
dt
dx
x
1
1
1
1
2
2
2
1
2
1
2
2
3
1
2





 

 
Substituindo na equação original, temos: 
uxaxaxa
dt
dx
a nn
n
n   10211  
Então, o sistema transformado passa a ser: 




















n
nnn
n
n
a
xaxaxau
dt
dx
x
dt
dx
x
dt
dx
x
dt
dx
12110
1
3
2
2
1
 

 
 
Em notação matricial: 
Teoria de Controle e suas aplicações 
 
225 
 u
a
x
x
x
x
x
x
a
a
a
a
a
a
a
a
a
a
a
a
dt
dx
dt
dx
dt
dx
dt
dx
dt
dx
dt
dx
n
n
n
n
n
n
n
n
n
n
nnn
n
n
n
















































































































































1
0
0
0
0
0
100000
010000
000000
000100
000010
1
2
3
2
1
123210
 
 
1
 
2
 
3
 
2
 
1









 
 
Podemos observar que a matriz A é formada por elementos que são iguais a 1 na 
posição imediatamente acima da diagonal principal, os elementos da última linha são 
iguais aos coeficientes da equação diferencial em sentido inverso dividido pelo primeiro 
coeficiente com o sinal trocado. Todos os demais elementos são nulos. 
15.8 Exemplos de aplicação dos métodos de integração 
Neste item apresentamos alguns exemplos de aplicação dos métodos de 
integração apresentados. 
15.8.1 Exemplo 1 
Sendo dado o seguinte circuito elétrico: 
 
Teoria de Controle e suas aplicações 
 
226 
a) Escrever o sistema de equações diferenciais de 1a ordem, utilizando apenas as 
variáveis de estado. 
b) Calcular os valores estimados para o passo de integração, o tempo máximo de 
simulação a serem utilizados na solução pelo método de Euler Modificado, o número 
de passos e o coeficiente de amortecimento. 
c) Indicar como resolver o sistema pelo método de Euler Modificado. 
d) Traçar os gráficos das variáveis de estado em função do tempo, destacando os seus 
respectivos valores iniciais e finais. 
Condições de contorno: chave aberta por um longo tempo. 
Solução: 
Encontrando a função 
dt
d Li 
L1C
L
0 iRe 
i 
 L e 
dt
d 
 L
iRee
 
i L1C0L 
dt
d
 (15.8.1.1) 
Encontrando a função 
dt
d Ce 
2
CC
L R
ee 
 Ci 
dt
d 
CR
eiR
 
e 
2
CL2C 
dt
d
 (15.8.1.2) 
Com a formulação em espaço de estados: 
uxy
uxx
 D C
 BA 


 
Teoria de Controle e suas aplicações 
 
227 


 

  
 
u
xy
y
y
u
x
x
dt
d
dt
d
0
C
L
2
1
0
C
L
2
1
C
L
e
D
0
0
e
i
C
10
01
e
B
0
L
1
e
i
A
CR
1
C
1
 L
1-
L
R
e 
i 













































































   
 
LC
1
LCR
R
e 
L
1i 
CR
1
CR
1
C
1
L
1-
L
R
CR
1e 
L
1-i 
i
1
2
1
CL
2
2
1
2
C
L
L







dt
d
dt
d
dt
d
dt
d
 
dt
d
dt
d CL
2
L
2
1 e 
L
1i 
CR
1
i 
LC
1
LCR
R






 
dt
d
dt
d CL
2
L
2
1 e i 
CR
L
i 
C
1
CR
R






 (15.8.1.3) 
Derivando a equação (15.8.1.1): 
0
i 
R
e i 
 L L1
C
2
L
2

dt
d
dt
d
dt
d (15.8.1.4) 
Substituindo a equação (15.8.1.3) na equação (15.8.1.4): 
0
i 
R
i 
CR
L
i 
C
1
CR
Ri 
 L L1
L
2
L
2
1
2
L
2







dt
d
dt
d
dt
d 
0
i 
CRR
i 
Li )R(R
i 
 LCR L21
L
L212
L
2
2  dt
d
dt
d
dt
d 
Teoria de Controle e suas aplicações 
 
228 
0i )R(R
i 
 C)RR(L
i 
 LCR L21
L
212
L
2
2  dt
d
dt
d
 (15.8.1.5) 
A equação é homogênea de 2ª ordem, com o seguinte polinômio característico: 
0 RR C)RR(L LCR 2121
2
2  rr 
As raízes do polinômio são: 
LC2R
)RLC(RR4C)RR(LC)RR(L
2
212
2
2121
21

,r 
 12,7268 2,833321 jjr ,  
Onde: 
Constante de tempo: s,
,
 35290
83332
11
T 

 
Período de oscilação: s,, 49370726812
2


 
Passo de integração: s,
,
n
,
h 01760
20
35290)min(




 (n=20 para Euler Modificado) 
Tempo máximo de simulação: s,tmax 76471)max(5  
Número de passos de integração: 100
h
t
n maxp 
Coeficiente de amortecimento: %,
,,
,
 7321
72681283332
83332
2222







 
 
 
 
 
 
 
Teoria de Controle e suas aplicações 
 
229 
Algoritmo para a solução pelo Método de Euler Modificado: 
 






















































































































 Façado Fim
 
enquanto Faça do Fim
2
BA
A 
2
BA
A 
1faça então 
e 
e e 
 OU 
i 
i i 
 E0 Se 
)B(A
2
ee 
)B(A
2
ii 
CR
)A(e)A(iR
 B
 L
)A(iR)A(ee
 B
Aee 
Aii 
0 
1 enquantoFaça 
 
1 
CR
eiR
A 
 L
iRee
A 
 2 1Faça 
0 t,0e ,0i :iniciais Valores
1
22
2
11
1
1
C
C
1
C
1
L
L
1
L
22C
1
C
11L
1
L
2
2C1L2
2
1L12C0
1
2CC
1LL
2
CL2
2
L1C0
1
11
C
1
L
htt
conv,toltolconv
h
h
hh
hh
h
h
conv
conv
conv
h
t
,,,p
pp
p
antp
p
antp
pp
pp
pp
pp
pant
pant
pp
pp
max
 
 
 
Teoria de Controle e suas aplicações 
 
230 
Gráficos: 
 
Valores iniciais: 
iL = 0 A 
eC = 0 V 
 
Valores finais: 
V 98040
051
50
e
RR
R
e
A 019610
051
1
RR
e
i
0
21
2
C
21
0
L
,
,










 
 
Diagrama em blocos no domínio da frequência: 
 
Tomando a transformada de Laplace da equação (15.8.1.1): 
Teoria de Controle e suas aplicações 
 
231 
  C0L1L1C0LL1C0L eei LRiReei L L
iRee
 
i 


 ss
dt
d
 
Então: 
 L R
ee
i
1
C0
L s

 . 
Tomando a transformada de Laplace da equação (15.8.1.2): 
  L2C2CL2C2
2
CL2C iRe CR1ei R e CR
CR
eiR
 
e 


 ss
dt
d
 
Então: 
C R 1
iR
 e
2
L2
C s
 . 
Com as expressões de IL e eC no domínio da frequência, podemos desenhar o 
seguinte diagrama em blocos: 
 
Matlab 
clear all 
close all 
tol=1.e-6 
itmax=100 
e0=1 
R1=1 
R2=50 
L=0.2 
C=0.03 
AA=[-R1/L -1/L; 
 1/C -1/(R2*C)] 
BB=[1/L; 
 0] 
CC=[1 0; 
Teoria de Controle e suas aplicações 
 
232 
 0 1] 
DD=[0; 
 0] 
lamb=eig(AA) 
a=R2*L*C 
b=L+R1*R2*C 
c=R1+R2 
r1=(-b+sqrt(b^2-4*a*c))/2/a 
r2=(-b-sqrt(b^2-4*a*c))/2/a 
sigma=real(lamb(1)) 
omega=imag(lamb(1)) 
T=-1/sigma 
Tal=2*pi/omega 
h=min(T,Tal)/20 
tmax=5*T 
np=tmax/h 
zeta=-sigma/sqrt(sigma^2+omega^2)*100 
step(AA,BB,CC,DD) 
t(1)=0; 
il(1)=0; 
ec(1)=0; 
for p=1:tmax/h 
 A1=(e0-ec(p)-R1*il(p))/L; 
 A2=(R2*il(p)-ec(p))/R2/C; 
 ily=il(p); 
 ecy=ec(p); 
 it=0; 
 erro=1; 
 while (erro==1 && it < itmax) 
 it=it+1; 
 B1=(e0-(ec(p)+h*A2)-R1*(il(p)+h*A1))/L; 
 B2=(R2*(il(p)+h*A1)-(ec(p)+h*A2))/R2/C; 
 ilx=il(p)+h/2*(A1+B1); 
 ecx=ec(p)+h/2*(A2+B2); 
 if ( abs((ilx-ily)/ilx) >tol || abs((ecx-ecy)/ecx) >tol ) 
Teoria de Controle e suas aplicações 
 
233 
 ily=ilx; 
 ecy=ecx; 
 A1=(A1+B1)/2; 
 A2=(A2+B2)/2; 
 else 
 erro=0; 
 end 
 end 
 il(p+1)=ilx; 
 ec(p+1)=ecx;t(p+1)=t(p)+h; 
end 
figure 
subplot(2,1,1) 
plot(t,il,'b-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'il (A)','FontName','Times','FontSize',12 ) 
grid 
axis( [ 0 1.8 -0.2 0.3 ] ) 
subplot(2,1,2) 
plot(t,ec,'b-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'ec (V)','FontName','Times','FontSize',12 ) 
grid 
axis( [ 0 1.8 0 1.5 ] ) 
15.8.2 Exemplo 2 
Sendo dado o seguinte circuito elétrico: 
 
Teoria de Controle e suas aplicações 
 
234 
 
a) Escrever o sistema de equações diferenciais de 1a ordem, utilizando apenas as 
variáveis de estado. 
b) Calcular os valores estimados para o passo de integração, o tempo máximo de 
simulação a serem utilizados na solução pelo método de Runge-Kutta de 4ª ordem, o 
número de passos e o coeficiente de amortecimento. 
c) Indicar como resolver o sistema pelo método de Runge-Kutta de 4ª ordem. 
d) Traçar os gráficos das variáveis de estado em função do tempo, destacando os seus 
respectivos valores iniciais e finais. 
Condições de contorno: chave aberta por um longo tempo. 
Encontrando a função 
dt
d Li : 





 









 
 C
L
2
1
L1
L
2
C
L
L1
L
0 e 
i 
 L
R
R
iR
i 
 L
R
e 
i 
 L
i R
i 
 Le
dt
d
dt
ddt
d
dt
d 
C
2
1
L
2
21L
2
21
C
2
1
L1
L
2
1
0 eR
R
i
R
RRi 
 L
R
RR
e
R
R
iR
i 
 L
R
R
1 e 




 







dt
d
dt
d 
 L)R(R
eRiRReR
 
i 
21
C1L2102L



dt
d
 (15.8.2.1) 
Encontrando a função 
dt
d Ce 
1
L
0
L
C
R
i 
 Le
i
e 
 C dt
d
dt
d 
 
dt
d
dt
d L
0L1
C
1
i 
 LeiR
e 
 CR  
Substituindo a expressão de 
dt
d Li 
Teoria de Controle e suas aplicações 
 
235 
)R(R
eRiRReR
e
 L)R(R
eRiRReR
 LeiR
e 
 CR
21
C1L2102
0
21
C1L2102
0L1
C
1 





dt
d 
C1L2102210L211
C
211 eRiRReR)R(Rei )R(RR
e 
 C )R(RR 
dt
d 
C1L2110L21L
2
1
C
211 eRiRRReiRRiR
e 
 C )R(RR 
dt
d 
C110L
2
1
C
211 eRReiR
e 
 C )R(RR 
dt
d 
C )R(R
iRee
 
e 
21
L1C0C



dt
d (15.8.2.2) 
Em formulação de espaço de estados: 
uxy
uxx
 D C
 B A


 

 

  
 
u
xy
y
y
u
x
x
dt
d
dt
d
0
C
L
2
1
0
21
21
2
C
L
2121
1
21
1
21
21
C
L
e
D
0
0
e
i
C
10
01
e
B
C )R(R
1
 L)R(R
R
e
i
A
C )R(R
1
C )R(R
R
 L)R(R
R
 L)R(R
RR
e 
i 









































































     
 
21
211
C
1
L
2
21
21
2
1
C
21
1L
21
2
21
2
1
2
21
21
C
21
1L
21
2121
1
21
1
21
21
21
C
21
1L
L
RR
)R(RR
e 
CR
i 
L
LC)R(R
RRR
e 
 LC)R(R
CRi 
 LC)R(R
L
LC)R(R
R
LC)R(R
RR
e 
 L)R(R
Ri 
C )R(R
1
C )R(R
1
C )R(R
R
 L)R(R
R
 L)R(R
RR
C )R(R
1e 
 L)R(R
Ri 
i






























dt
d
dt
d
dt
d
dt
d
dt
d
dt
d
dt
d
dt
d
 
dt
d
dt
d CL
1
L
e 
C
i 
R
L
i  
Teoria de Controle e suas aplicações 
 
236 
0
e 
CR
i 
 Li R C1
L
L1  dt
d
dt
d (15.8.2.3) 
Derivando a equação (1): 
dt
d
dt
d
dt
d C
1
L
212
L
2
21
e 
R
i 
RR
i 
 L)R(R  
dt
d
dt
d
dt
d L
22
L
2
1
21C i R
i 
 L
R
RRe 


 (15.8.2.4) 
Substituindo a equação (4) na equação (3): 
0
i 
R
i 
 L
R
RR
CR
i 
 Li R L22
L
2
1
21
1
L
L1 








dt
d
dt
d
dt
d 
0
i 
CRR 
i 
 LC)R(R
i 
LiR L212
L
2
21
L
L1  dt
d
dt
d
dt
d 
0iR
i 
C)RR(L
i 
 LC)R(R L1
L
212
L
2
21  dt
d
dt
d 
0i 
RR
Ri 
RR
CRRLi 
LC L
21
1L
21
21
2
L
2






dt
d
dt
d (15.8.2.5) 
 
Equação homogênea de segunda ordem com o seguinte polinômio característico: 
0
RR
R
RR
CRRL
 LC
21
1
21
212 




 rr 
 
As raízes do polinômio são: 
2LC
RR
LCR4
RR
CRRL
RR
CRRL
21
1
2
21
21
21
21
21













,r 
 jjr ,  9,3314 2,857121 
 
Teoria de Controle e suas aplicações 
 
237 
Constante de tempo: s,
,
 35000
85712
11


 
Período de oscilação: s,
,
 67330
33149
22



 
Passo de integração: s,
,
n
),min(
h 035000
10
35000



 (n=10 para Runge-Kutta de 4ª 
ordem) 
Tempo máximo de simulação: s,)max(tmax 75015  
Número de passos de integração: 50
h
t
n maxp 
Coeficiente de amortecimento: %,
,,
,
 27729
3314985712
85712
2222







 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Teoria de Controle e suas aplicações 
 
238 
Método de Runge-Kutta de 4ª ordem: 
 
 
 






















































 




 







 




 







 




 







 




 














h
h
h
hh
hh
hh
hh
hh
hh
htp
pp
pp
pp
pp
pp
pp
pp
pp
pp
pp
pp
max
t t
)D2C2B(A
6
ee 
)D2C2B(A
6
ii 
C )R(R
)C(iR)C(ee
 D
 L)R(R
CeR)C(iRReR
 D
C )R(R
2
B
iR
2
B
ee
C 
 L)R(R
2
B
eR
2
B
iRReR
C 
C )R(R
2
A
iR
2
A
ee
 B
 L)R(R
2
A
eR
2
A
iRReR
 B
C RR
iRee
A 
 L)R(R
eRiRReR
A 
/ 1, 
0t
0e
0i
1
2222C
1
C
1111L
1
L
21
1L12C0
2
21
2C11L2102
1
21
1
L1
2
C0
2
21
2
C1
1
L2102
1
21
1
L1
2
C0
2
21
2
C1
1
L2102
1
21
L1C0
2
21
C1L2102
1
1
1
C
1
L
 
 
 
 
 
Teoria de Controle e suas aplicações 
 
239 
Gráficos: 
 
 
Valores iniciais: 
iL = 0 A 
eC = 0 V 
Valores finais: 
iL = e0 / R1 = 0,05 A 
eC = 0 V 
 
Matlab 
 
clear all 
close all 
Teoria de Controle e suas aplicações 
 
240 
e0=1 
R1=20 
R2=1 
L=0.2 
C=0.05 
AA=[-R1*R2/(R1+R2)/L R1/(R1+R2)/L; 
 -R1/(R1+R2)/C -1/(R1+R2)/C] 
BB=[R2/(R1+R2)/L; 
 1/(R1+R2)/C] 
CC=[1 0; 
 0 1] 
DD=[0; 
 0] 
lamb=eig(AA) 
a=L*C 
b=(R1*R2*C+L)/(R1+R2) 
c=R1/(R1+R2) 
r1=(-b+sqrt(b^2-4*a*c))/2/a 
r2=(-b-sqrt(b^2-4*a*c))/2/a 
sigma=real(lamb(1)) 
omega=imag(lamb(1)) 
T=-1/sigma 
Tal=2*pi/omega 
h=min(T,Tal)/10 
tmax=5*T 
np=tmax/h 
zeta=-sigma/sqrt(sigma^2+omega^2)*100 
step(AA,BB,CC,DD) 
t(1)=0; 
il(1)=0; 
ec(1)=0; 
for p=1:tmax/h 
 A1=(R2*e0-R1*R2*il(p)+R1*ec(p))/((R1+R2)*L); 
 A2=(e0-ec(p)-R1*il(p))/(R1+R2)/C; 
 B1=(R2*e0-R1*R2*(il(p)+h*A1/2)+R1*(ec(p)+h*A2/2))/((R1+R2)*L); 
Teoria de Controle e suas aplicações 
 
241 
 B2=(e0-(ec(p)+h*A2/2)-R1*(il(p)+h*A1/2))/(R1+R2)/C; 
 C1=(R2*e0-R1*R2*(il(p)+h*B1/2)+R1*(ec(p)+h*B2/2))/((R1+R2)*L); 
 C2=(e0-(ec(p)+h*B2/2)-R1*(il(p)+h*B1/2))/(R1+R2)/C; 
 D1=(R2*e0-R1*R2*(il(p)+h*C1)+R1*(ec(p)+h*C2))/((R1+R2)*L); 
 D2=(e0-(ec(p)+h*C2)-R1*(il(p)+h*C1))/(R1+R2)/C; 
 il(p+1)=il(p)+h/6*(A1+2*B1+2*C1+D1); 
 ec(p+1)=ec(p)+h/6*(A2+2*B2+2*C2+D2); 
 t(p+1)=t(p)+h; 
end 
figure 
subplot(2,1,1) 
plot(t,il,'b-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'il (A)','FontName','Times','FontSize',12 ) 
grid 
axis( [ 0 1.8 0 0.08 ] ) 
subplot(2,1,2) 
plot(t,ec,'b-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'ec (V)','FontName','Times','FontSize',12 ) 
grid 
axis( [ 0 1.8 -0.04 0.08 ] ) 
 
15.8.3 Exemplo 3 
Sendo dado o seguinte circuito elétrico: 
Teoria de Controle e suas aplicações 
 
242 
 
 
a) Escrever o sistema de equações diferenciais de 1a ordem, utilizando apenas as 
variáveis de estado. 
b) Calcular os valores estimados para o passo de integração, o tempo máximo de 
simulação a serem utilizados na solução pelo método de Runge-Kutta de 4ª ordem, o 
número de passos e o coeficiente de amortecimento. 
c) Indicar como resolver o sistema pelo método de Runge-Kutta de 4ª ordem.d) Traçar os gráficos das variáveis de estado em função do tempo, destacando os seus 
respectivos valores iniciais e finais. 
 
Condições de contorno: chave aberta por um longo tempo tempo e capacitor C2 
descarregado antes do fechamento da chave. 
 
Encontrando a função 
dt
d Li 
) cos(EiR 
i 
 Le mL1
L
c1   tdt
d 
 L
iRe E
 
i L1C1mL 
)tcos(
dt
d 
 (15.8.3.1) 
Teoria de Controle e suas aplicações 
 
243 
Encontrando a função 
dt
d C1e 
2
C2C1C1
1L R
eee 
 Ci


dt
d
 
12
L2C1C2C1
CR
iRee
 
e 

dt
d
 (15.8.3.2) 
Encontrando a função 
dt
d C2e 
C2C1
C2
22 ee
e 
CR 
dt
d
 
22
C2C1C2
CR
ee
 
e 

dt
d
 (15.8.3.3) 
 
Em formulação de espaço de estados: 
uxy
uxx
 D C
 B A


 
 

 

 
 
u
xy
y
y
y
u
x
x
dt
d
dt
d
dt
d
0
C2
C1
L
2
2
1
0
C2
C1
L
2222
12121
1
C2
C1
L
e
D
0
0
0
e
e
i
C
100
010
001
e
B
0
0
L
1
e
e
i
A
CR
1
CR
1
0
CR
1
CR
1
C
1
0
L
1
L
R
e 
e 
i 









































































































   
 
 
 
 
Teoria de Controle e suas aplicações 
 
244 
Método de Runge-Kutta de 4ª ordem: 
 
Parâmetros: 
000010 0000010 010 100 50
0 
3
2
138000 2 60
2121 ,C,C,LR,R
Eff m

 
 
 
Valores iniciais: 
0
0
1
1
tan
2
cos
1
1
 
tancos
1
1
2
2
1
2
11
1
11
1
1
2
1
2
1
1
11
1
































 
































 




t
e
C
LRC
R
C
L
E
e
C
LR
R
C
L
E
i
C
m
C
m
L










 
Teoria de Controle e suas aplicações 
 
245 
     
     
   










































































 




 






 




 




 






 




 


 






 




 






 




 




 






 




 


 












ht t
DCBA
h
e e
DCBA
h
e e
DCBA
h
ii
CR
hCehCe
 D
CR
hBiRhBehBe
D
L
hCiRhCe)htE
D
CR
B
he
B
he
C
CR
B
hiR
B
he
B
he
C
L
B
hiR
B
he
h
tω E
C
CR
A
he
A
he
B
CR
A
hiR
A
he
A
he
B
L
A
hiR
A
he
h
tω E
B
CR
ee
A
CR
iRee
A
iRetE
A
htp
pp
p
C
p
C
p
C
p
C
p
L
p
L
p
C
p
C
p
L
p
C
p
C
p
L
p
C
p
m
p
C
p
C
p
L
p
C
p
C
p
L
p
C
p
m
p
C
p
C
p
L
p
C
p
C
p
L
p
C
p
m
p
C
p
C
p
L
p
C
p
C
p
L
p
C
p
m
1
33332
1
2
22221
1
1
1111
1
22
3221
3
12
122132
2
1121
1
22
3
2
2
1
3
12
1
2
2
1
3
2
2
1
1
2
1
1
22
3
2
2
1
3
12
1
2
2
1
3
2
2
1
1
2
1
1
22
21
3
12
212
2
11
1
max
)22(
6
)22(
6
)22(
6
 
 
 
( cos
 
22
 
222
 
 
22
)
2
(cos
 
22
 
222
 
 
 
22
)
2
(cos
 
 
 
 L
) cos(
 
/ 1, 




 
Teoria de Controle e suas aplicações 
 
246 
 
Gráficos: 
Preto: eC1 e iL 
cinza: eC2 e iC2 
 
 
Alterando-se o valor de C1 para o valor de C2, pode-se fazer uma simulação de operação 
de fechamento em back-to-back, e os gráficos ficam: 
 
Preto: eC1 e iL 
cinza: eC2 e iC2 
 
Teoria de Controle e suas aplicações 
 
247 
 
 
Matlab 
 
clear all 
close all 
 
e0=138000*sqrt(2)/sqrt(3) 
f=60 
w=2*pi*f 
fi=0 
R1=0.5 
R2=100 
L=0.01 
C1=0.000001 
C2=0.00001 
A=[-R1/L -1/L 0; 
Teoria de Controle e suas aplicações 
 
248 
 1/C1 -1/R2/C1 1/R2/C1; 
 0 1/R2/C2 -1/R2/C2] 
B=[1/L; 
 0 
 0] 
C=[1 0 0; 
 0 1 0; 
 0 0 1] 
D=[0; 
 0; 
 0] 
 
y=[ 1 1/(R2*C1)+1/(R2*C2)+R1/L R1/(R2*L*C1)+R1/(R2*L*C2)+1/(L*C1) 
1/(R2*L*C1*C2) ] 
roots(y) 
 
lamb=eig(A) 
sigma1=real(lamb(1)) 
sigma2=real(lamb(3)) 
omega=imag(lamb(1)) 
freq=omega/2/pi 
T1=-1/sigma1 
T2=-1/sigma2 
Tal=2*pi/omega 
h=min(T1,T2); 
h=min(h,Tal)/10 
Tmax=5*max(T1,T2) 
Tf=1/f 
Tmax=max(Tmax,2*Tf) 
np=Tmax/h 
zeta1=-sigma1/sqrt(sigma1^2+omega^2)*100 
zeta2=100 
step(A,B,C,D) 
t(1)=0; 
il(1)=e0*cos(fi-atan((w*L-1/w/C1)/R1))/sqrt(R1*R1+(w*L-1/w/C1)^2); 
Teoria de Controle e suas aplicações 
 
249 
ec1(1)=e0*cos(fi-pi/2-atan((w*L-1/w/C1)/R1))/sqrt(R1*R1+(w*L-1/w/C1)^2)/w/C1; 
ec2(1)=0; 
ic2(1)=0; 
for p=1:Tmax/h 
 a1=(e0*cos(w*t(p)+fi)-R1*il(p)-ec1(p))/L; 
 a2=(ec2(p)-ec1(p)+R2*il(p))/R2/C1; 
 a3=(ec1(p)-ec2(p))/R2/C2; 
 b1=(e0*cos(w*(t(p)+h/2)+fi)-R1*(il(p)+h*a1/2)-(ec1(p)+h*a2/2))/L; 
 b2=((ec2(p)+h*a3/2)-(ec1(p)+h*a2/2)+R2*(il(p)+h*a1/2))/R2/C1; 
 b3=((ec1(p)+h*a2/2)-(ec2(p)+h*a3/2))/R2/C2; 
 c1=(e0*cos(w*(t(p)+h/2)+fi)-R1*(il(p)+h*b1/2)-(ec1(p)+h*b2/2))/L; 
 c2=((ec2(p)+h*b3/2)-(ec1(p)+h*b2/2)+R2*(il(p)+h*b1/2))/R2/C1; 
 c3=((ec1(p)+h*b2/2)-(ec2(p)+h*b3/2))/R2/C2; 
 d1=(e0*cos(w*(t(p)+h)+fi)-R1*(il(p)+h*c1)-(ec1(p)+h*c2))/L; 
 d2=((ec2(p)+h*c3)-(ec1(p)+h*c2)+R2*(il(p)+h*c1))/R2/C1; 
 d3=((ec1(p)+h*c2)-(ec2(p)+h*c3))/R2/C2; 
 il(p+1)=il(p)+h/6*(a1+2*b1+2*c1+d1); 
 ec1(p+1)=ec1(p)+h/6*(a2+2*b2+2*c2+d2); 
 ec2(p+1)=ec2(p)+h/6*(a3+2*b3+2*c3+d3); 
 ic2(p+1)=(ec1(p+1)-ec2(p+1))/R2; 
 t(p+1)=t(p)+h; 
end 
ilmin=min(il) 
ilmax=max(il) 
ec1min=min(ec1) 
ec1max=max(ec1) 
ec2min=min(ec2) 
ec2max=max(ec2) 
ic2min=min(ic2) 
ic2max=max(ic2) 
ecmin=min(ec1min,ec2min); 
ecmax=max(ec1max,ec2max); 
imin=min(ilmin,ic2min); 
imax=max(ilmax,ic2max); 
figure 
Teoria de Controle e suas aplicações 
 
250 
subplot(2,1,1) 
hold on 
plot(t,ec1,'b-','LineWidth',2) 
plot(t,ec2,'r-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'ec1, ec2 (V)','FontName','Times','FontSize',12 ) 
box 
grid 
hold off 
axis( [ 0 Tmax ecmin ecmax ] ) 
subplot(2,1,2) 
hold on 
plot(t,il,'b-','LineWidth',2) 
plot(t,ic2,'r-','LineWidth',2) 
xlabel( 'Tempo (s)','FontName','Times','FontSize',12 ) 
ylabel( 'il, ic2 (A)','FontName','Times','FontSize',12 ) 
box 
grid 
axis( [ 0 Tmax imin imax ] ) 
hold off 
 
15.8.4 Exemplo 4 
Resolver o sistema abaixo pelo método trapezoidal implícito: 
 
Teoria de Controle e suas aplicações 
 
251 
 
 
Ml
lKMgsen
dt
d
lKMgsenFPsen
dt
d
Ml
dt
dv
MMaF
lKKvF
lv
aR
a
 
 )(
 
 









 
 










 
M
K
sen
l
g
dt
d
dt
d
 
 
0 
 
2
2
2
2




sen
l
g
dt
d
M
K
dt
d
M
K
sen
l
g
dt
d
dt
d
 
 
com: 
g = 9,81 m/s2 
Teoria de Controle e suas aplicações 
 
252 
l = 1 m 
K = 0,5 kg/s 
M = 2 kg 
 
Considerar  = 0 e  = 60º como condições iniciais. 
 
Fazer os gráficos de  versus t,  versus t e  versus . 
 
Para pequenas amplitudes de   sen , então: 
0 
2
2
 
l
g
dt
d
M
K
dt
d
 
Cuja equação característica é 0 2 
l
g
M
K  . 
Então os autovalores são 
2
4
 
2
l
g
M
K
M
K






 , que geralmente são 
complexos conjugados, e com isso temos uma constante de tempo T e um período de 
oscilação , dados por: 
s 8
2

K
M
T e s 0082
4
4
2
,
M
K
l
g








 
Se for desprezado o atrito: s 00622 ,
g
l
  
Teoria de Controle e suas aplicações 
 
253 
 
Figura 15.9 – Posição angular 
 
 
Figura 15.10 – Velocidade angular 
 
Teoria de Controle e suas aplicações 
 
254 
 
Figura 15.11 – Trajetória  ×  
 
 
Figura 15.12 – Força de tração na barra