ENGEST 014 cee14

Prévia do material em texto

ISSN 1413-9928 
(versão impressa) 
 
CADERNOS DE 
ENGENHARIA DE ESTRUTURAS 
Universidade de São Paulo 
Escola de Engenharia de São Carlos 
Departamento de Engenharia de Estruturas 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
______________________________________________________________________ 
 
Uma família de algoritmos hermitianos para a integração 
direta das equações de dinâmica das estruturas 
 
 
Heitor Miranda Bottura 
José Elias Laier 
 
______________________________________________________________________ 
 
 
 
 
 
 
 
 
Número 14 
São Carlos, 1999 
 
 
UMA FAMÍLIA DE ALGORITMOS HERMITIANOS PARA A INTEGRAÇÃO 
DIRETA DAS EQUAÇÕES DE DINÂMICA DAS ESTRUTURAS 
 
Heitor Miranda Bottura1 & José Elias Laier2 
 
 
 
RESUMO 
 
 Apresenta-se aqui uma família de algoritmos de passo simples, com ordem de 
precisão local qualquer e aniquilamento assintótico para a análise dinâmica de 
estruturas. São utilizadas expressões hermitianas para as relações em diferenças 
envolvidas na representação das equações que descrevem o problema. Explicitam-se os 
membros da família, com precisão desde a primeira até a oitava ordem, que apresentam 
estabilidade incondicional, efetuando-se sua análise espectral. Um exemplo com dois 
graus de liberdade é resolvido, com o objetivo de mostrar o efeito do aniquilamento 
assintótico sobre a representação do movimento. 
 
 
 
1. INTRODUÇÃO 
 
 A solução do problema de análise dinâmica de estruturas semidiscretizadas no 
espaço sob excitação qualquer variável com o tempo passa pela solução numérica de um 
sistema de equações diferenciais ordinárias de segunda ordem, que expressa as 
condições de equilíbrio segundo o sistema de coordenadas utilizado. 
 A solução analítica fechada de tais equações é um caminho que logo se esgota, por 
limitar bastante a natureza da solicitação. Outras formas de abordar o problema são em 
geral necessárias. Algumas delas, como o emprego das transformadas de Fourier e da 
integral de Duhamel também não são suficientemente flexíveis, ou demandam trabalho 
numérico excessivo, para satisfazer o analista na maioria dos casos, a exemplo dos 
problemas não-lineares. 
 A solução mais geral é a que substitui as relações diferenciais do problema por 
expressões finitas, gerando um sistema de equações lineares simultâneas onde os dados 
de entrada configuram o movimento em instantes já determinados e as incógnitas são 
representadas pelo deslocamento e suas derivadas num instante posterior, por 
determinar. 
 Dada a diversidade de solicitações de interesse na análise estrutural, dificilmente 
se chegará a um método ideal para um caso genérico, sendo conveniente estudar-se o 
maior número possível de suas características, para melhor definir sua viabilidade e 
possibilidades de emprego. 
 O primeiro grupo, ou família de métodos, a merecer destaque é o dos processos 
derivados da aplicação de expressões em diferenças: especificamente a utilização de 
expressões em diferenças finitas centrais e o método apresentado por HOUBOLT 
 
1 Prof. Assistente do Departamento de Engenharia Civil, UNESP, Bauru. E-mail: 
heitor@bauru.unesp.br 
2 Prof. Titular do Departamento de Engenharia de Estruturas, EESC-USP. E-mail: jelaier@sc.usp.br 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
2
(1950), uma bem conhecida variante dessa família. HULBERT (1991) apresenta uma 
generalização deste último, que é aperfeiçoada por HULBERT (1994) e HULBERT & 
CHUNG (1994b). Sua importância maior é a de esgotar as possibilidades dos métodos 
de passo simples de ordens inferiores, assim entendidos os que apresentam erro local de 
até segunda ordem. 
 O processo de NEWMARK (1959) introduz a idéia de parâmetros livres cuja 
variação acarreta modificações nos coeficientes das expressões que permitem 
intensificar características consideradas de utilidade. A idéia é retomada por WILSON 
(1968) com o seu método Wilson-θ. Posteriormente , outros autores apresentaram 
processos com parâmetros livres, como BRUSA & NIGRO (1980), WOOD, BOSSAK 
& ZIENKIEWICZ (1980), BAZZI & ANDERHEGGEN (1982) com o que 
denominaram família ρ e KATONA & ZIENKIEWICZ (1985), com o seu método Beta-
m. 
 Uma alternativa à de escrever expressões em diferenças é a formulação por 
resíduos ponderados, introduzida por ZIENKIEWICZ (1977), que discretiza o domínio 
do tempo similarmente à forma utilizada para o tratamento de domínios elásticos no 
espaço por elementos finitos. Este caminho tem seqüência em ZIENKIEWICZ, WOOD 
& TAYLOR (1980) e HOFF & PAHL (1988) entre outros. 
 O uso de mínimos quadrados é outra abordagem que se prestou à montagem de 
algoritmos, como no caso da família de KUJAWSKI & GALLAGHER (1989) . 
 O tratamento simultâneo do domínio espaço-tempo por elementos finitos é a 
proposta de HUGHES & HULBERT (1988) e este caminho é retomado por HULBERT 
(1992 e 1994), aproveitando o método descontínuo de Galerkin (TDG). 
 A apresentação da formulação hermitiana aplicada à solução de equações 
diferenciais ordinárias de primeira ordem é descrita em COLLATZ (1968) e está 
presente nas análises efetuadas nos trabalhos de MAKINSON (1968) e NØRSETT 
(1973). Posteriormente, expressões hermitianas de primeira e terceira ordem são 
utilizadas por ARGYRIS, VAZ & WILLAM (1977) na solução de problemas de 
escoamento viscoso e de viscoelasticidade, viscoplasticidade e difusão transiente. Uma 
extensão da técnica à equação de equilíbrio dinâmico já é feita por ARGYRIS & al. 
(1973), e essa idéia é retomada posteriormente por ARGYRIS & al (1979), inclusive 
com a presença do aniquilamento assintótico. Em 1991, ARGYRIS & MLEJNEK 
voltam ao tema, agora sem obter dissipação numérica. 
 Um algoritmo hermitiano cúbico é gerado por LAIER (1993), com parâmetros 
livres de ajuste, e pode ser comparado favoravelmente ao de KUJAWSKI & 
GALLAGHER (1989). Podem-se citar ainda os estudos publicados por LAIER (1995). 
 
 
2. OBJETIVOS 
 
 O objetivo do presente trabalho é a obtenção de algoritmos para a integração 
numérica da equação diferencial de equilíbrio, empregada na representação do 
movimento de estruturas sob solicitações dependentes do tempo. 
 Partiu-se da formulação hermitiana para as expressões em diferenças, a serem 
utilizadas no lugar das relações diferenciais presentes no problema. 
 Estabeleceram-se requisitos a serem atendidos pelos métodos resultantes, com 
vistas a que fossem competitivos comparativamente aos já existentes. 
 A ordem de precisão local deveria ser maior do que dois, de modo a configurar o 
que a literatura chama de “método de ordem superior”, uma vez que a bibiliografia 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
3
citada na introdução deste trabalho aponta para o esgotamento dos métodos de ordens 
um e dois, também chamados de métodos de ordens inferiores. 
 Propriedades tradicionalmente reputadas como importantes para o comportamento 
do algoritmo deveriam estar presentes. Consideraram-se nessa relação a consistência, 
convergência e estabilidade, essa última com ênfase nos casos em que independe do 
tamanho do passo adotado para o estudo do movimento, isto é, quando é incondicional. 
 Dado que é importante que um algoritmo se preste à análise nãolinear, métodos 
multipassos foram descartados, intentando-se chegar apenas a processos de passo 
simples. 
 O aniquilamento assintótico, que permite filtrar os modos superiores aos 
correspondentes a uma certa “freqüência de corte”, eliminando sua influência, foi outra 
exigência fixada.Não se estabeleceu a necessidade da presença de parâmetros livres de ajuste, por 
vezes utilizados para dosar ou fazer aparecer características desejadas nos algoritmos. 
 A forma escolhida de estudar-se as propriedades dos métodos resultantes foi, 
como tradicionalmente feito na literatura, a análise espectral efetuada na aplicação a 
sistemas com um único grau de liberdade, sob vibrações livres sem amortecimento. 
 Um exemplo de sistema com mais de um grau de liberdade deveria ser 
apresentado, a fim de destacar a ação do aniquilamento assintótico. 
 
 
3. A EQUAÇÃO DE EQUILÍBRIO 
 
 Uma estrutura cuja descrição do movimento exija o conhecimento do 
deslocamento em mais de uma coordenada de um sistema qualquer escolhido diz-se de 
múltiplos graus de liberdade, e a equação matricial de equilíbrio segundo essas 
coordenadas tem a forma: 
 
 [ ]{ } [ ] [ ] ( ){ }tPxMxCxK =++
•••
 (1) 
 
onde {x} é o vetor dos deslocamentos nas coordenadas adotadas, um ponto indica uma 
derivação com relação ao tempo e dois pontos indicam a segunda derivada de {x}, isto 
é, a aceleração. A matriz de massa é escrita [M], o amortecimento é representado por 
[C] e o vetor das solicitações, por {P(t)}. A matriz de rigidez [K] do sistema é aqui 
assumida constante. 
 Embora se vislumbre a aplicação dos métodos obtidos à solução de problemas 
nãolineares, para se chegar a conclusões definitivas numa análise espectral trata-se aqui, 
numa primeira abordagem, do caso linear, como é tradicional. 
 A solução de (1) fornece os deslocamentos no sistema de coordenadas escolhido, 
arranjados no vetor {x} para o instante considerado, de modo que se tenha {x}T = {x1, 
x2, ..xn}. O índice n corresponde aqui portanto ao número de graus de liberdade do 
sistema. 
 Ao se estudarem as propriedades de um algoritmo de integração, é vantajoso 
considerar sistemas com um único grau de liberdade (SDOF) e a literatura discute a 
validade da extrapolação das conclusões obtidas sobre o método. Para associar 
simplicidade de tratamento matemático e generalidade às conclusões, aborda-se 
inicialmente o caso de vibrações livres. Assim, (1) passa a escrever-se: 
 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
4
 0xmxckx =++ ••• (2) 
 
na nomenclatura descrita anteriormente. É ainda usual estabelecerem-se as relações: 
 
 
c
ncn c
c= ; m2c ; 
m
k γω==ω 
 
de modo que agora o equilíbrio se escreva na forma: 
 
 0xx2x 2nn =ω+γω+
•••
 (3) 
 
mais sintética, onde ωn representa a freqüência natural de vibração do movimento, cc é o 
amortecimento crítico e γ é a chamada taxa de amortecimento do movimento. 
 A expressão geral de (3), segundo a representação utilizada por COLLATZ 
(1968), é: 
 
 0x,x,x,tF =

 ••• (4) 
 
válida para qualquer t pertencente ao domínio de F. Assim, tomado um ti, de 
deslocamento e derivadas conhecidas, como referência, o equilíbrio para um instante ∆t 
à sua frente, isto é, distante de um passo, pode ser escrito então: 
 
 0x,x,x,tF 1i1i1i =

∆ +••+•+ (5) 
 
ou, alternativamente: 
 
 

∆= +••+•+ 1i1i1i x,x,tfx (6) 
 
expressão que relaciona as três incógnitas de interesse para avançar no movimento. 
 
 
4. ALGORITMOS DE INTEGRAÇÃO 
 
 As derivadas em (6) podem ser substituídas por relações finitas envolvendo 
valores já determinados que indiquem a tendência do movimento, e que tenham 
geralmente aspecto: 
 
 
( )1rg1i1i1iiii1i1i1i 1tR,...x,x,x,x,x,x,x,x,tgx +−••−•−•••+••++• ∆+∆= (7) 
 
( )1rh1i1i1iiii1i1i1i 2tR+,...x,x,x,x,x,x,x,x,thx +−••−•−•••+•++•• ∆∆= (8) 
 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
5
onde a escolha das expressões que combinam os termos conhecidos para a obtenção dos 
deslocamentos e derivadas incógnitos define o método empregado. Pode-se mostrar que 
é sempre possível escrever (7) e (8) de ordem de aproximação r não-negativa, isto é, tais 
expressões finitas diferem das derivadas que representam por um polinômio em ∆t de 
coeficientes não-nulos apenas para expoentes superiores a r, onde r1 e r2 são as ordens 
de aproximação das representações respectivamente de velocidade e aceleração. 
 O conjunto de (6), (7) e (8), omitidos os resíduos, reproduz um sistema de 
equações lineares nas incógnitas 1i1i1i x,x,x +
••
+
•
+ cuja solução permite avançar um passo 
no tempo. 
 Conforme o arranjo das relações e seus coeficientes, é possível resolver o sistema 
sem a inversão de matrizes, isto é, não se trata de um sistema simultâneo de equações; 
basta uma multiplicação de matrizes para resolvê-lo sem necessidade de inverter 
nenhuma delas. O ganho computacional é evidente e o método direto é chamado de 
explícito. Métodos em que isso não é possível são chamados de implícitos e essa 
diferenciação tem conseqüências práticas nas características resultantes no processo. 
 O método explícito mais tradicional é o das diferenças centradas e o atrativo da 
eficiência computacional nos processos dessa natureza, principalmente nos problemas 
com não-linearidades, tem estimulado a busca de novos algoritmos dessa classe, como 
atestam os trabalhos de CHUNG (1994) e HULBERT & CHUNG (1994b) e PEZESHK 
(1995) entre outros, bem como a apresentação de versões "explícitas" que certos 
algoritmos permitem quando se atendem certas condições, como é feito em BAZZI & 
ANDERHEGGEN (1982). 
 A maior desvantagem desses algoritmos é a impossibilidade de apresentarem 
estabilidade incondicional como KRIEG (1973) mostra, exigindo controle do tamanho 
do passo para evitar a deterioração da resposta. 
 Quando é suficiente utilizar informações relativas apenas ao passo anterior, o 
algoritmo resultante é dito de passo simples; caso contrário, trata-se de um 
procedimento de passo múltiplo. O fato de ser imediata a alteração do passo nos 
algoritmos de passo simples é uma vantagem importante e a solução de problemas não-
lineares praticamente exige métodos dessa natureza. 
 Para requintar os conhecimentos sobre o problema podem-se utilizar relações que 
envolvam derivadas superiores. Para a explicitação das novas incógnitas em ti+1, torna-
se necessário levantar a indeterminação gerada, o que pode ser feito derivando-se (6) 
com relação ao tempo tantas vezes quantas preciso, pois as expressões resultantes são 
independentes. Assim, o sistema de equações é complementado com as relações: 
 
 
( )
( ) ( ) 

∆==


∆==
+
••
+++
••
+
•
++
••
+
•
1i
4
1i
3
1i2
2
1i
1i
3
1i1i1i
p,x,x,tf
dt
fdx
p,x,x,tf
dt
dfx
 (9) 
 
e assim por diante. Evidentemente, há que pensar também numa representação finita 
para as derivadas da excitação nos movimentos forçados. 
 A obtenção das informações no instante ti utiliza o mesmo caminho, isto é, 
derivações sucessivas da equação do equilíbrio,com a mesma observação. 
 Tais algoritmos cobram um preço pela possibilidade de obtenção de melhores 
resultados; para estruturas com n deslocamentos generalizados incógnitos, ou graus de 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
6
liberdade resultantes da semidiscretização, geram sistemas de ordem 2n de equações 
lineares, o dobro comparativamente aos processos que envolvem derivadas até segunda 
ordem apenas. 
 Os métodos assim descritos são chamados "de ordens superiores" e sua utilização 
foi inicialmente preterida até que a aplicação de técnicas especiais de solução de 
sistemas aumentaram sua competitividade. 
 
 
5. CARACTERÍSTICAS DE ALGORITMOS DE INTEGRAÇÃO 
 
5.1. Auto-iniciação 
 
 Há processos cuja aplicação não é possível desde o início do movimento, isto é, 
não bastam as condições iniciais e a escolha de um passo para a obtenção do 
deslocamento e derivadas no instante t1 diretamente a partir da utilização das relações 
que constituem o método. É o caso dos algoritmos multipassos, como se pode observar 
nas equações (7) e (8), tomando-se nelas i = 0. 
 Nesses casos, para viabilizar a implementação do processo, aplica-se um 
procedimento de partida, para tantos instantes iniciais quantos necessários para a 
completa montagem das relações que caracterizam o algoritmo, de modo que as 
propriedades resultantes da solução numérica podem estar contaminadas por esse 
procedimento. 
 
5.2. Estabilidade 
 
 O conceito de estabilidade numérica que interessa utilizar é aquele associado à 
garantia de que as aproximações para os deslocamentos e suas derivadas produzidas 
pelo método não cresçam indefinidamente, para o tempo tendendo ao infinito. 
 Seja um algoritmo definido por um conjunto de relações como (6), (7) e (8), 
omitindo-se o resíduo e rearranjado para assumir o aspecto seguinte: 
 
 
[ ] [ ]






∆
∆−=






∆
∆
••
•
+
••
+
•
+
i
2
i
i
1i
2
1i
1i
xt
xt
x
N
xt
xt
x
M
 (10) 
 
ou, fazendo: - [M]-1[N] = [A] e representando mais sinteticamente: 
 
 { } [ ]{ }i1i xAx =+ (11) 
 
onde se indica a obtenção do deslocamento e suas derivadas no instante ti+1 como 
resultado de uma transformação aplicada aos respectivos valores no instante ti. Tomado 
o instante ti+n, (11) é escrita: 
 { } [ ] { }inni xAx =+ (12) 
 
onde se observa, fazendo n crescer indefinidamente, que a estabilidade é governada pelo 
comportamento da chamada matriz [A] de amplificação do método. 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
7
 Quanto à possível influência das condições iniciais na estabilidade, a dúvida é 
dirimida a partir da observação da relação (11) escrita para i = 0. De fato, para qualquer 
{x0} tal propriedade apresenta-se dependente apenas da matriz [A]. 
 Seja a decomposição espectral de [A], dada por: 
 
 [ ] [ ][ ] [ ] 1nn PPA −Λ= 
 
onde [P] é a matriz dos autovetores ortonormalizados e [Λ], uma matriz diagonal 
contendo os autovalores λi, i=1,2,...m de [A]. Seja ρ(A) o raio espectral de [A], cuja 
definição é dada pela expressão: 
 
 ( ) m,...2,1i,máxA i =λ=ρ (13) 
 
Então, de acordo com BATHE (1976), [Λ]n é limitado para n tendendo ao infinito se, e 
somente se, for atendida a condição: 
 
 ( ) 0,1A ≤ρ (14) 
 
que define o critério de estabilidade empregado para classificar os algoritmos. 
 Quando, para atender (14) é necessário tomar um passo menor do que um certo 
valor, o processo é dito condicionalmente estável. Passos maiores têm, como 
conseqüência, menor esforço computacional. É comum a utilização de intervalos 
maiores que os necessários para garantir estabilidade nos modos superiores dos 
algoritmos condicionalmente estáveis. É, portanto, interessante dispor-se de estabilidade 
incondicional. 
 
5.3. Consistência e precisão 
 
 Seja um algoritmo, que possa ser representado na forma: 
 
 { } [ ]{ } ( ){ }2rii1i tRxAx ++ ∆+= (15) 
 
onde r+2, correspondente à ordem de precisão nos deslocamentos, é o menor expoente 
de ∆t no vetor {Ri} dos resíduos. 
 Segundo HILBER & HUGHES (1978), consistência é a propriedade apresentada 
pelos métodos em que r é positivo. 
 
5.4. Convergência 
 
 Se quanto menor o passo, mais próxima estiver a resposta obtida do valor correto, 
o algoritmo é chamado de convergente, como definem BAZZI & ANDERHEGGEN 
(1982). Invocando o teorema de Lax-Richtmayer, HILBER* apud autores acima afirma 
 
* HILBER, H. M. Analysis and design of numerical integration methods in structural dynamics. Report 
n. EERC 76-29, Berkeley, Earthquake Engineering Research Center, University of California, 1976 
apud BAZZI, G. ; ANDERHEGGEN, E. The ro-family of algorithms for time-step integration 
with improved numerical dissipation. Earthqu. Eng. Struct. Dyn., v. 10, p. 537-550, 1982. 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
8
que algoritmos consistentes e estáveis no limite para o passo tendendo a zero são 
convergentes. 
 Essas condições acima enunciadas podem ser expressas por meio do par de 
desigualdades: 
 
 
( )


>
≤ρ→∆ω
0r 
0,1Alim 0t
 (16) 
 
e são suficientes inclusive para vibrações forçadas, desde que a excitação p(t) seja uma 
função analítica, quando o teorema de existência para equações a coeficientes analíticos 
garante a convergência da solução de (3) em séries de potências. 
 
5.5. Aniquilamento 
 
 O erro resultante das operações numéricas efetuadas na solução dos modos 
superiores de um movimento tende a ser maior do que o presente nos primeiros modos. 
Ademais, sua parcela contribuinte na representação do deslocamento e suas derivadas é, 
em geral, pouco significativo, embora isto possa não ser verdadeiro nos problemas 
envolvendo impacto. 
 Visando anular a influência desses modos superiores, alguns algoritmos 
introduzem uma dissipação numérica seletiva, que os atinge após os instantes iniciais. 
 Essa propriedade estará presente nos métodos que respeitarem o limite dado por: 
 
 ( ) 0Alim t =ρ∞→∆ω (17) 
 
e pode ser vista como a inclusão de um amortecimento numérico encarregado de 
aniquilar os modos superiores. 
 
5.6. Overshoot 
 
 Alguns algoritmos estáveis superestimam os deslocamentos (e/ou derivadas) nos 
instantes iniciais do movimento. A explicação para isso é que o comportamento nos 
primeiros instantes é governado não pelo raio espectral mas pela norma da matriz de 
amplificação do algoritmo. 
 A primeira observação dessa propriedade indesejável foi feita por GOUDREAU & 
TAYLOR (1972) ao analisar o algoritmo WILSON-θ. Foi-lhe dada a denominação de 
"overshoot",numa referência a seu efeito direto de superestimar os primeiros 
deslocamentos e suas derivadas. 
 
 
 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
9
5.7. Erro global e erro local 
 
 Fazendo-se referência apenas aos erros de discretização ou trucamento, o erro 
global cometido no deslocamento é dado pela diferença entre o valor fornecido pelo 
método empregado e o deslocamento teórico exato. Via de regra, não é possível o 
cálculo do valor acima pelo desconhecimento deste último. 
 Tomem-se os valores iii x,x,x
•••
 do deslocamento e derivadas aproximados 
fornecidos pelo método para ti, como condições iniciais de um movimento, e 
represente-se por xi+1 a solução exata correspondente para o instante ti+1. A 
expressão para o erro local escreve-se neste caso: 
 ( ) 1i1i1i xxxE +++ −=l (18) 
 
isto é, decorre do avanço em um simples passo do movimento. Seu valor tampouco 
pode ser calculado no caso geral, porém sua ordem é definida pela ordem de precisão r 
do algoritmo utilizado. 
 Uma diferença importante entre os dois conceitos enunciados pode ser entendida 
com o exemplo de um método que apresente pequeno erro local, porém sistemático, 
como um acréscimo no valor de uma grandeza a cada passo. O acúmulo desse pequeno 
erro local resultará, após um certo número de passos, num erro global considerável. 
 Infelizmente, o mais difícil de ser avaliado ou estimado é justamente o erro global; 
a previsão da qualidade do resultado obtido é feita lançando-se mão de extrapolações de 
análises de movimentos teoricamente conhecidos ou da influência da ordem do erro 
local sobre o erro global. 
 Para os métodos de passo simples, se o erro local é de ordem r+1, HAIRER & 
WANNER* mostram que o erro global é uma combinação deste mais os erros de 
transporte, e é de ordem r. 
 
5.8. Dissipação e dispersão 
 
 Apesar das dificuldades enunciadas de medir os erros de discretização há valores 
tradicionalmente referidos pela bibliografia para comparar algoritmos de integração. 
 Tomam-se geralmente vibrações livres e nãoamortecidas, para condições iniciais 
correspondentes a um deslocamento da posição de repouso x(t=0) = x0, sem velocidade 
inicial. 
 Este movimento tem como vantagens principais ser periódico, facilitando a 
medida de erros numéricos nas freqüências ou períodos resultantes, além de ter 
amplitude teórica constante, sendo portanto imediata a quantificação do erro no seu 
cálculo. Além disso a ausência de excitação evita associar-se conclusões à natureza da 
solicitação. 
 A equação característica de (10), associada a este movimento, tem o aspecto: 
 
 [ ] [ ]( ) 0IAdet =λ−− (19) 
 
cujos coeficientes reais são os invariantes da matriz de amplificação. A possibilidade de 
serem todas reais pode ser verificada a partir do critério de Routh-Hurwitz, citado em 
 
* HAIRER, E.; NØRSETT, S. P. ; WANNER, G. Solving ordinary differential equations 1: Nonstiff 
problems. Berlin, Springer, 1987. 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
10
KRIEG (1973). É de se esperar que as condições necessárias não sejam atendidas, dadas 
as características do movimento que representa de acordo com as hipóteses assumidas. 
Resta, portanto, a possibilidade de um par de raízes complexas conjugadas λ1 e λ2, e 
uma terceira raíz real, λ3, chamada "espúria". 
 Para os algoritmos convergentes, de muito maior interesse prático, existe uma 
constante positiva ~Ω tal que, se tomado um passo de modo que Ω = ω∆t pertença ao 
intervalo (0,~Ω) pode-se afirmar a validade da relação: 
 
 
0,12,13 ≤λ<λ (20) 
 
indicando-se que as raízes complexas λ1 e λ2 são mais importantes que λ3 na 
representação do movimento. 
 Adotando-se a forma polar para representar as duas raízes principais, sendo ρ seu 
módulo, utilizando o índice num para referir-se aos valores resultantes da aplicação do 
método numérico, e considerando θ = ω∆t, bem como a nomenclatura: 
 
 θ
ρ−=γ lnnum
 ; t
num ∆
θ=ω
 (21) 
 
pode-se escrever a expressão do deslocamento alternativamente, na forma: 
 
 ( ) ( )[ ]tnsenCtncosCex num2num1tnn numnum ∆ω+∆ω= ∆ωγ− (22) 
 
porque a literatura mostra que, mesmo quando λ3 não é significativamente inferior à 
unidade, pouca relevância tem na composição da resposta numérica. De qualquer forma, 
basta ter em mente que a influência de λ3 desaparece em algoritmos convergentes para 
o limite do passo tendendo a zero, o que justifica a omissão do termo que lhe 
corresponde de ora em diante. 
 Comparando-se (22) com a solução exata de (3) sem o termo do amortecimento, 
dada por: 
 
 tcosAtsenAx n2n1 ω+ω= (23) 
 
vê-se a introdução de um fator formalmente similar ao existente em problemas com 
amortecimento físico. Isto é decorrente da dissipação numérica produzida pelo método, 
dada por (22). Algoritmos sem aniquilamento obviamente não têm dissipação. 
 O amortecimento numérico γnum é uma medida proporcional ao erro global 
presente. Entretanto, sua presença controlada é desejável, por aniquilar os modos 
superiores. 
 Note-se o compromisso entre o erro global e aniquilamento produzido; por tal 
razão, comparações de métodos sempre remetem ao cotejo da dissipação neles presente. 
 Observa-se, também, uma discordância no período fornecido como resposta, pois: 
 
 n
2T ω
π=
 ; num
num
2T ω
π=
 
 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
11
e pode-se utilizar o erro relativo do período, dado por: 
 
 
1
T
TT
num
nnum −ω
ω=−=τ
 (24) 
 
como uma medida da distorção da freqüência, chamada de dispersão. 
 
 
6. A FORMULAÇÃO HERMITIANA 
 
 A expressão hermitiana para diferenças, na forma de passo simples, escreve-se: 
 
 ( ) ( ) ( )[ ] ( )∑
=
+ =∆++∆
n
0j
r
i
j
1ij
j
ij
j 0tRxbxat (25) 
e a substituição dos valores assumidos no ponto i+1 por séries de Taylor centradas em i 
levadas em (25) resulta em: 
 
 
( ) ( )
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )( )
( ) ( )( ) ( ) 0tR..xtbxta
..txxtbxta
..x
!2n
t..txxtbxta
...x
!1n
t..x
2
txtxtbxta
..x
!n
t..x
6
tx
2
txtxbxa
r
i
n
i
n
n
n
i
n
n
n
i
1n
i
1n
1n
1n
i
1n
1n
n
i
2n
3
ii
2
2i
2
2
n
i
1n
3
i
2
ii1i1
n
i
n
3
i
3
i
2
ii0i0
=∆++∆+∆
++∆+∆+∆+
+


 +−
∆++∆+∆+∆+
+


 +−
∆++∆+∆+∆+∆+
+


 +∆++∆+∆+∆++
−−
−
−−
−
−••••
−••••
•••
 (26) 
 
e agrupando-se os coeficientes das derivadas de mesma ordem em ∆t, obtêm-se as 
identidades: 
 
 
( )
( )
( ) ( ) 0xtba..!1n
b
!n
b
.
0xtbab
2
b
0xtbab
0xban
i
n
nn
10
i
2
221
0
i110
i00
=∆


 +++−+
=∆

 +++
=∆++
=+
••
•
 (27) 
 
em número total de n+1. Essas relações geram um sistema de equações homogêneas nos 
2(n+1) coeficientes ak e bk. Pode-se tomar o conjunto nãonulo e arbitrar-se valor para 
um deles, o que fornece para o grau de indeterminação do sistema o valor n. 
 Este último corresponde ao número máximo p de parâmetros livres presentes em 
(25) para que essa expressão em diferenças se anule até os termos daquela ordem em ∆t. 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
12
 Se considerados nãonulos e a determinar, os coeficientes ak e bk, k = 0,1,..m que 
respeitem (27) até n valem as expressões: 
 
 p + (n + 1) = 2m + 1 ⇒ n = 2m - p (28) 
 
resultando ordem igual a n + 1 para o erro, ou seja, observando-se a expressão acima 
vale: 
 
 r = 2m - p + 1 . (29) 
 
 A expressão do equilíbrio para vibrações livres nãoamortecidas é escrita: 
 
 x xi i
••
+ ++ =1 2 1 0ω (30) 
 
para o instante seguinte ao de referência do movimento. Se tomados os operadores 
hermitianos para representar (7) e (8), mas utilizando-se as relações que (30) e suas 
derivadas permitem escrever, pode-se adotar a representação: 
 
 
( )
( ) 0tRx,x,x,x,tH
0tRx,x,x,x,tH
2r2
i1i1iii
2
n
2r1
i1i1iii
1
n
2
2
1
1
=∆+

∆
=∆+

∆
++
•
+
•
++
•
+
•
 (31) 
 
onde ∆t é o passo escolhido, ti é um instante de referência de deslocamento e derivadas 
já conhecidos e ainda utiliza-se a notação ti+1 = ti + ∆t como um instante onde essas 
grandezas são desconhecidas e a determinar. 
 A segunda parcela corresponde ao resíduo ou erro de truncamento local cometido 
no cálculo da expressão de diferenças que substitui as derivadas presentes. Sua ordem é 
r1 + 2 no primeiro algoritmo das (31), onde a ordem da maior derivada considerada em 
(25) para a sua montagem é dada por n1. 
 Se os algoritmos forem independentes, trata-se de um sistema de duas equações 
lineares simultâneas nas incógnitas xi+1 e xi
•
+1 , cuja solução permite o avanço de um 
passo no conhecimento do movimento. 
 Dessa forma, para um sistema com n graus de liberdade, a ordem do sistema 
resultante é 2n. Isto é característico dos métodos de ordens superiores. Visando-se 
gerar uma família de algoritmos que apresente aniquilamento assintótico, observa-se 
que para tanto deve-se obedecer a (17), que é atendida se: 
 
 2,1j,i 0alim ijt ==∞→∆ω (32) 
 
ou seja, se todos os elementos da matriz de amplificação tenderem a se anular com o 
aumento indefinido do passo. 
 Efetuando-se as operações necessárias, os elementos de [A] escrevem-se: 
 
 [ ] ( )
( ) ( )
( ) ( )


−−
−−
−−= 1221221111212111
2212122221121122
21122211 nmnmnmnm
nmnmnmnm
mmmm
1A (33) 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
13
 
com a expressão do resíduo dada analogamente por: 
 
 ( ){ } ( ) ( )( )



∆
∆


−
−
−=∆ +
+
+
2r2
2r1
1121
1222
21122211
2r
2
1
tR
tR
mm
mm
mmmm
1tR (34) 
 
também em função dos coeficientes dos deslocamentos e derivadas envolvidos. 
 Observando-se (10) verifica-se que os mij são funções dos coeficientes do 
deslocamento e suas derivadas no instante ti, enquanto os nij dependem de valores 
associados ao instante ti+1. 
 Se levada (30) a (26), pode-se também reescrever cada um dos operadores para 
vibrações livres nãoamortecidas, explicitando-se os primeiros termos da somatória 
como: 
 
 
0...xbtxatxtbtxtat
xbtxatxtbxtaxbxa
1i4
44
i4
44
1i3
22
i3
22
1i2
22
i2
22
1i1i11i0i0
≈+∆ω+∆ω+

∆∆ω−

∆∆ω−
+∆ω−∆ω−∆+∆++
++
••
++
••
+
 
 
ou agrupando-se os termos relativos ao instante de referência e ao seguinte, ∆t à sua 
frente: 
 
 
( ) ( )
( ) ( ) 0xt..btbx..btb
xt..atax..atata
1i3
22
11i2
22
0
i3
22
1i4
44
2
22
0
≈∆+∆ω−++∆ω−+
+∆+∆ω−+−∆ω+∆ω−
+
•
+
•
 (35) 
 
de modo que se este for o primeiro operador considerado, juntando-se o segundo vem: 
 
 0xtmxmxtnxn 1i121i11i12i11 ≈∆++∆+ +
•
+
•
 
 0xtmxmxtnxn 1i221i21i22i21 ≈∆++∆+ +
•
+
•
 (36) 
 
desde que se omitam as parcelas referentes ao resíduo. 
 Objetivando-se atender a (38) e considerando-se (39) bem como (42) e (43), é 
conveniente tomar a subfamília de operadores representados por: 
 
 ( ) ( ) ( ) ( ) ( ) 0tRxbtxatH n
0j
2r
i
m
0k
k
1ik
kj
ij
j
n,m =∆+∆+∆=∑ ∑
=
+
=
+ (37) 
isto é, tomar valores independentes de m e n que garantam a mij ordem em ∆t grande o 
suficiente quando comparada à de nij para conduzir ao desejado aniquilamento. 
 O número de coeficientes a determinar é m + n + 2; fixado um não-nulo e 
arbitrário para relacioná-lo aos demais, tem-se como número de equações necessárias m 
+ n + 1. 
 Este último número exprime também a ordem de anulação do erro de truncamento 
nos deslocamentos, isto é: 
 
 r + 2 = m + n + 1 ⇒ r = m + n - 1 (38) 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
14
 
sendo r a ordem de precisão local nos esforços e acelerações. 
 Pode-se sacrificar a precisão em favor de parâmetros de ajuste visando-se a 
alguma propriedade. Então passa a valer: r = m + n - 1 - p, sendo p o número de 
parâmetros. 
 Para a montagem de um algoritmo que respeite as condições acima indicadas 
como desejáveis, um critério para a escolha do segundo operador é tomar como 
expressão a derivada do primeiro com relação ao tempo. Com isso, os coeficientes são 
os mesmos e conseqüentemente é de esperar que as características sejam mantidas, tal 
como a ordem relativa (em ∆t) dos coeficientes de mij e nij, que são um indicativo de 
que o algoritmo é mais ou menos “conservador”. 
 Tomando-se, então, o operador: 
 
 ( ) ( ) ( ) ( ) ( ) 0tRxbtxatH 2r1m
0k
k
1ik
k
n
0j
j
ij
j1
n,m
1 =∆+∆+∆= +
=
+
=
∑∑ (39) 
considere-se o segundo, atendendo ao proposto acima, definido por: 
 
 
( ) ( ) ( ) ( ) ( ) ( )∑ ∑+
=
+
=
+
+++ =∆+∆+∆==
1n
1j
1m
1k
3r2k
1ik
kj
ij
j
1
n,m2
1n,1m 0tRxdtxctdt
Hd
H 1 (40) 
 
onde valem as relações decorrentes da derivação cj = aj-1 e dk = bk-1, resultando: 
 
 ( )( ) ( ) ( ) ( ) ( ) 0tRxbtxatH 3r21m
1k
k
1i1k
k
1n
1j
j
i1j
j2
1n,1m
1 =∆+∆+∆= +
+
=
+−
+
=
−++ ∑∑ . (41) 
 
 A substituição da equação do movimento, sempre para vibrações livres 
nãoamortecidas, em (41) e (39), à semelhança de (35) mas tomando-se θ = ω∆t como 
tradicionalmentefeito na literatura, fornece então a forma final do algoritmo dada por: 
 
 
[ ] [ ]
[ ] [ ]
[ ] [ ]
[ ] [ ] 0R...aaaxt...aaax 
...bbbxt...bbbxH
0R...aaaxt...aaax+ 
...bbbxt...bbbxH
2
i4
4
2
2
0i5
6
3
4
1
2
i
4
4
2
2
01i5
6
3
4
1
2
1i
2
1n,1m
1
i5
4
3
2
1i4
4
2
2
0i
5
4
3
2
11i4
4
2
2
01i
1
n,m
=+−θ+θ−∆++θ−θ+θ−+
+−θ+θ−∆++θ−θ+θ−=
=+−θ+θ−∆+−θ+θ−
+−θ+θ−∆+−θ+θ−=
•
+
•
+++
•
+
•
+
 (42) 
 
e de uma comparação com (36) pode-se escrever que valem para os elementos das 
matrizes envolvidas as expressões que seguem: 
 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
15
 
1122
12
2
21
7
6
5
4
3
2
112
6
6
4
4
2
2
011
1122
12
2
21
7
6
5
4
3
2
112
6
6
4
4
2
2
011
nn
nn
...aaaan
...aaaan
mm
mm
...bbbbm
...bbbbm
=
θ−=
+θ−θ+θ−=
+θ−θ+θ−=
=
θ−=
+θ−θ+θ−=
+θ−θ+θ−=
 (43) 
 
e considerando-se (10) para este caso particular, tem-se para o determinante de [M] o 
valor: 
 
 [ ] 2122211 mmMdet θ+= (44) 
 
completamente definido a partir dos bk. 
 
 
7. AS PROPRIEDADES ESPECTRAIS 
 
 Para os elementos da matriz de amplificação [A] valem as relações: 
 
 
[ ]
[ ]
12
2
21
11121211
12
22
1212
2
1111
11
aa
Mdet
nmnma
a
Mdet
nmnma
θ−=
−−=
=θ+−=
 (45) 
 
levando-se em conta (35) e sua derivada. 
 A equação característica de [A] assume aspecto: 
 
 ( ) 0aa 2122211 =θ+λ− 
 
resultando o par de raízes complexas conjugadas dado por: 
 
 1211 aia θ±=λ 
 
e consequentemente com raio espectral expresso por: 
 
 2
12
22
11
2
12
22
112
12
22
11 mm
nnaa θ+
θ+=θ+=λ=ρ 
 
resultado a que se pode chegar considerando-se as relações desenvolvidas 
anteriormente. 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
16
 É interessante notar que a bifurcação de raízes está na origem. De fato, apenas no 
caso limite de θ nulo ocorre uma raiz real dupla, o que elimina a preocupação com 
raízes ao se efetuar um estudo de desempenho do processo em questão. 
 Levando-se em conta (43) tem-se a expressão final para o raio espectral dessa 
família: 
 
 
( ) ( )
( ) ( )27654321226644220
2
7
6
5
4
3
2
1
22
6
6
4
4
2
2
0
..bbbb..bbbb
..aaaa..aaaa
+θ−θ+θ−θ++θ−θ+θ−
+θ−θ+θ−θ++θ−θ+θ−=ρ (46) 
 
que é um ponto de partida para o estudo das propriedades dos processos em pauta. 
 Algumas características são inerentes aos algoritmos de passo simples, como a 
auto-iniciação. 
 As derivadas são aproximadas a partir de sua expansão em séries de Taylor, do 
que resultam métodos sempre consistentes. Sua expressão é semelhante a (15), com r1 e 
r2 positivos. 
 A convergência está garantida uma vez que sempre se terá |a0| = |b0|, respeitando-
se assim a condição imposta por (16). 
 Da forma como se montou o algoritmo, de (46) verifica-se a existência do 
aniquilamento assintótico desde que se tome m maior que n, pois isto resultará em 
denominador com ordem em ∆t superior à do numerador, atendendo-se também a (17); 
a ordem do erro alcançada é dada pela expressão (38). 
 Um estudo da estabilidade pode ser feito de forma imediata, confrontando-se (14) 
com a expressão (46) resultante do algoritmo. 
 A comparação de desempenho com outros métodos é tradicionalmente feita 
medindo-se a dissipação numérica introduzida, dada por (21) bem como a distorção 
resultante, assim denominado o erro relativo no período, explicitado pela relação (24). 
 Apresentam-se, a seguir, os primeiros membros dessa família e ilustram-se suas 
propriedades e características. 
 
 
8. OS ALGORITMOS DA FAMÍLIA 
 
 A determinação do primeiro algoritmo assim obtido sem a inclusão de parâmetros 
livres de ajuste é feita a partir da consideração de (38) tomado r unitário e impondo-se 
aniquilamento assintótico, o que formece: 
 
 m + n = 2 e m > n 
 
havendo assim a única possibilidade de m = 2 e n = 0, correspondendo a: 
 
 
( )
( ) ( ) ...xt
6
1tR
0tRx
2
txtxxH
3
i
331
i
31
i1i
2
1i1ii
1
0,2
+∆−=∆
=∆+∆−∆+−= +••+•+
 
 
combinado com a sua derivada relativamente ao tempo, de aspecto: 
 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
17
 
( ) ( )
( ) ( ) ...xt
6
1tR
0tRx
2
txtxtxtH
4
i
442
i
42
i
3
1i
3
1i
2
1ii
2
1,3
+∆−=∆
=∆+∆−∆+∆−∆= ++
••
+
••
 
 
reproduzindo assim o método de primeira ordem apresentado por LAIER (1995). 
 Para atender à estabilidade incondicional, deve-se respeitar a: 
 
 ( ) ( )[ ] θ∀≤ρ⇒≤ρ ;0,1A0,1A 2 . 
 
 De (46), substituindo-se os valores de ak e bk encontrados, resulta: 
 
 ( ) ( )
( )
00,1
1
2
11
01 4
22
2
2
222
2 ≥θ⇒≤
θ+

 θ+−
θ+=ρ 
 
verdadeiro para qualquer θ, tratando-se portanto de um algoritmo incondicionalmente 
estável. A matriz de amplificação é definida pela aplicação das (44) e (45), com os 
valores: 
 a0 = 1, b0 = -1, b1 = 1 e b2 = - 1/2 
 
representando os coeficientes não-nulos presentes na combinação linear representativa 
das expressões em diferenças que reproduzem as derivadas no tempo presentes no 
problema. Para o algoritmo de segunda ordem, as condições impostas são: 
 
 m + n = 3 e m > n 
 
o que leva a duas possibilidades. Na primeira delas, para n unitário e m = 2, o operador 
resultante é: 
 
 
( )
( ) ( ) ...xt
12
1tR
0tRxtxt4xt2x6x6H
4
i
441
i
41
i1i
2
1ii1ii
1
1,2
+∆=∆
=∆+∆+∆−∆−+−= +••+••+
 
 
já utilizado por LAIER (1995) mas não conjuntamente com o seu derivado, de aspecto: 
 
( ) ( )
( ) ( ) ...xt
12
1tR
0tRxtxt4xt2xt6xt6H
5
i
552
i
52
i
3
1i
3
1i
2
i
2
1ii
2
2,3
+∆=∆
=∆+∆+∆−∆−∆+∆−= ++
••••
+
••
. 
 
cuja verificação de estabilidade em (46) escrita com os ak e bk encontrados, produz: 
 
 ( ) ( )( ) ( )( ) 00,1246
26 4
222
222
2 ≥θ⇒≤−θ+θ−
−θ+−=ρ 
 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
18
essa exigência é atendida incondicionalmente para qualquer passo adotado. 
 A evolução das demais propriedades com o tamanho do passo é mostrada nas 
figuras 1, 3 e 5. 
 Sua matriz de amplificação define-se pelas expressões (44) e (45) tomando-se: 
 
 a0 = -6, b0 = 6, a1 = -2, b1 = -4, b2 = 1 
 
que são os valores encontrados para os coeficientes da expressão em diferenças. 
 O outro operador possível é o correspondente a n nulo e m = 3. Encontra-se a 
relação: 
 
 
( ) ( )
( ) ( ) ...xt
4
1tR
0tRxtxt3xt6x6x6H
4
i
441
i
41
i
3
1i
3
1i
2
1i1ii
1
0,3
+∆=∆
=∆+∆+∆−∆+−= ++
••
+
•
+
 
 
a ser utilizado em conjunto com 2 1,4H para que o método resultante pertença à família 
em questão. Entretanto,ao se verificar a estabilidade, nota-se a exigência de passos 
correspondentes a θ2 ≥ 3 para a presença dessa propriedade. Não se tenciona aqui 
abordar processos condicionalmente estáveis, razão pela qual se deixa de lado esse 
algoritmo. 
 Não se descarta, porém, a possibilidade de outras vantagens suplantarem esse 
"defeito" para a utilização em problemas específicos, o que poderia ser comprovado 
levando-se a efeito um estudo mais profundo da tentativa de implementação ao caso 
particular de interesse do analista. 
 Outro comentário pertinente refere-se ao fato de que tal operador não apresenta 
obstáculos de per si para seu emprego. Basta abandonar a exigência de pertencer à 
família em foco para que haja inúmeras possibilidades de combinações com outros 
operadores. Aliás, seu potencial "genético" é bastante bom, dadas as condições que se 
procurou atender por ocasião da escolha dos coeficientes considerados na expressão de 
diferenças. 
 Outra possibilidade de uso de tais operadores é a introdução, a partir de suas 
expressões básicas, de têrmos cujos coeficientes sejam parâmetros de ajuste que 
permitam restaurar a estabilidade incondicional, embora sacrificando-se a ordem de 
precisão. 
 Aplicando-se o procedimento descrito acima, podem-se determinar por inspeção, 
os membros incondicionalmente estáveis de ordens superiores da família. A tabela 1 
apresenta os valores dos coeficientes correspondentes aos métodos de primeira a oitava 
ordem. A última linha fornece o valor do coeficiente do termo de menor ordem do erro 
local. 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
19
 
Tabela 1 - Coeficientes dos termos dos algoritmos e do erro local correspondente 
 
ORDEM 1 2 3 4 5 6 7 8 
a0 1 -6 24 60 360 840 6720 15120 
b0 -1 6 -24 -60 -360 -840 -6720 -15120 
a1 -2 6 24 120 360 2520 6720 
b1 1 -4 18 36 240 480 4200 8400 
a2 3 12 60 360 1260 
b2 -1/2 1 -6 -9 -72 -120 -1200 -2100 
a3 4 20 120 
b3 1 1 12 16 200 300 
a4 5 
b4 -1 -1 -20 -25 
a5 
b5 1 1 
C.R. -1/6 1/4 1/20 -1/24 -1/210 1/1680 1/3024 1/30240 
 
 Os algoritmos descritos têm lei de formação igual e geral. Nada há a indicar que 
não se possam obter membros de ordens superiores às aqui apresentadas, e com as 
mesmas características, adotando-se o mesmo critério para descartar os "indesejáveis". 
Entretanto, o objetivo central, ao se explicitarem os primeiros métodos resultantes, é 
apenas a geração dos gráficos das principais propriedades. 
 
 
0
0,2
0,4
0,6
0,8
1
1,2
0,1 1 10 100 1000
WDT
RO1
RO2
RO3
RO4
 
Fig. 1 - Raio espectral dos algoritmos hermitianos de primeira a quarta ordem 
 
 
0
0,2
0,4
0,6
0,8
1
0,1 1 10 100 1000
WDT
R
O
RO5
R06
RO7
RO8
 
 
Fig. 2 - Raio espectral x passo dos algoritmos hermitianos de quinta a oitava ordem 
 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
20
0
2
4
6
8
10
0,00 0,40 0,80 1,20 1,60 2,00 2,40 2,80
WDT
G
AM
A*
10
0% GAMA1%
GAMA2%
GAMA3%
GAMA4%
 
Fig. 3 - Amortecimento numérico dos algoritmos hermitianos de primeira a quarta 
ordem 
 
 
0
2
4
6
8
10
0,00 0,80 1,60 2,40 3,20 4,00 4,80
WDT
GAMA5%
GAMA6%
GAMA7%
GAMA8%
 
Fig. 4 - Amortecimento numérico dos algoritmos hermitianos de quinta a oitava ordem 
 
 
 
0
0,1
0,2
0,3
0,4
0,5
0,00001 0,40001 0,80001 1,20001 1,60001 2,00001 2,40001 2,80001 3,20001 3,60001 4,00001
WDT
ALONG1
ALONG2
ALONG3
ALONG4
 
Fig. 5 - Dispersão x passo nos algoritmos hermitianos de primeira a quarta ordem 
 
 
 
0
0,1
0,2
0,3
0,4
0,5
0,0 0,4 0,8 1,2 1,6 2,0 2,4 2,8 3,2 3,6 4,0
WDT
ALONG5
ALONG6
ALONG7
ALONG8
 
Fig. 6 - Dispersão x passo nos algoritmos hermitianos de quinta a oitava ordem 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
21
 
9. UM EXEMPLO NUMÉRICO 
 
A figura 7 mostra o esquema estrutural de um pórtico de dois andares cujas vigas têm 
rigidez igual e unitária. Não há amortecimento físico no problema e as condições 
iniciais correspondem a deslocamento e velocidade nulos para ambos os graus de 
liberdade, numerados conforme indicado na figura. A massa do sistema é representada 
por uma unidade concentrada no nível de cada andar. 
 Este esquema é utilizado em um exemplo apresentado por WARBURTON (1976) 
porém sob outra solicitação. Nas condições descritas, as matrizes de massa e rigidez 
bem como a excitação para instantes a partir do início do movimento, saindo do 
repouso, são dadas respectivamente por: 
 
 
 
 
 
 
 
 
 
 
 
 
 
 Fig. 7 - Esquema estrutural para o segundo exemplo 
 
 


−
−=

=
11
12
]K[ ; 
10
01
]M[ 



−= 4,8
0,10
}P{ 
 
e efetuando-se a decomposição segundo os modos de vibração, tem-se: 
 
 s1664,10T ; 
8506,0
52157,0
}{ ; 382,0
2
53
11
2
1 ≅


=φ≅−=ω 
 
 s8832,3T ; 
5257,0
8506,0
}{ ; 618,2
2
53
22
2
2 ≅


−=φ≅+=ω . 
levando às soluções dadas por: 
 
 
).1t6180,1(cos9358,4y
)1t6181,0(cos9425,4y
2
1
−=
−=
 
 
 Quanto à solução estática, isto é, desprezando-se o efeito das acelerações, resta 
considerarem-se as relações: 
 
 [ ]{ } [ ] { }




−
−=



⇒Φ=Ω
94,4
94,4
y
y
)t(Py
ST,2
ST,1T2 
 
X1
X2
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
22
vetor que fornece os deslocamentos da estrutura. Em relação ao sistema proposto, vale: 
 
 { } [ ]{ }yx Φ= 
 
de forma que, após o aniquilamento do segundo modo, se pode escrever: 
 
 { }






 −=



=
ST,2
1
AN,2
AN,1
AN y
y
5257,08506,0
8506,05257,0
x
x
x 
 
expressão que fornece os deslocamentos quando se considera a contribuição do 
primeiro modo e se aniquilou completamente a influência do segundo. Pode-se dizer 
que esta é a resposta nos deslocamentos após a “filtragem” completa do segundo modo 
de vibração. Num problema de múltiplos graus de liberdade, raciocínio análogo é 
aplicável, relativamente aos modos inferiores e superiores, assim denominados a partir 
de sua posição relativamente à freqüência de corte ou “filtragem” pretendida. A 
experiência e sensibilidade do analista e natureza do problema são fatores que envolvem 
a determinação de uma adequada freqüência de corte de um movimento em um método 
direto que apresente aniquilamento assintótico. 
 Os deslocamentos correspondentes ao primeiro modo são dados por: 
 
 { } [ ]



Φ=



=
0
y
x
x
x 1
1M,2
1M,1
1M 
e, portanto, a presença do aniquilamento pode ser medida através de: 
 
 { } { } { } { } { }( ) { }Rxxxxx ST1MNUMANNUM =+−=− 
 
onde {xNUM} contém os deslocamentos calculados numericamente pelo algoritmo, e 
{R} é um vetor nas diferenças indesejadas entre os deslocamentos fornecidos pelo 
método e os pretendidos. Há duas contribuições em {R}: a do modo inferior amortecido 
prematuramente e a do modo superior, ainda eventualmente presente. 
 Um bom algoritmo deve ser capaz de apresentar {R} com forte tendência a se 
anular, mesmo quando as freqüências cuja influência se quer separar estiverem 
relativamente próximas. No caso em pauta, ω1 ≈ 0,6181 s-1 e ω2 ≈ 1,6180 s-1. A tabela 2 
mostra o resíduo parax1, resultante da aplicação dos algoritmos de terceira, quarta e 
quinta ordens da família apresentada, para um passo ∆t = T1 /4. 
 Observa-se, para a quarta ordem, que o resíduo praticamente se anula após vinte 
passos, o que significa fiel representação do modo inferior e aniquilamento 
praticamente total do superior. 
 A tabela mostra ainda que, para o algoritmo de quinta ordem, o aniquilamento 
tarda mais a efetuar-se. Isso pode ser compensado vantajosamente utilizando-se passos 
maiores, a fim de acelerar a presença do fenômeno. Os gráficos de raio espectral versus 
tamanho do passo apresentados neste trabalho mostram que o “corte” das freqüências 
superiores é cada vez mais abrupto, à medida que se tomam membros da família de 
ordem mais elevada. 
 O uso de ordens mais baixas pode levar à indesejável situação de, no caso de 
freqüências próximas, não se aniquilarem as superiores sem que o amortecimento 
numérico comece a atuar sobre os modos inferiores, no intervalo de observação do 
movimento. 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
23
 A coluna correspondente ao resíduo do algoritmo de terceira ordem mostra isso. 
Assim, intervalos de observação maiores estão associados à necessidade de métodos de 
ordens mais elevadas, dosando-se o aniquilamento a partir do tamanho de passo 
adequado. 
 O comportamento dos valores citados está representado nas figuras 6.8 a 6.10. 
 
Tabela 2 - Deslocamentos x1 e resíduo (R), para ∆t = T1 / 4 
 
N x1 M1 x1 x1 H3 R3 x1 H4 R4 x1 H5 R5
0 0 0 0 -4,202 0 -4,202 0 -4,202
1 -2,598 3,971 2,938 1,334 4,064 2,46 3,579 1,975
2 -5,197 0,523 -1,013 -0,018 -1,065 -0,07 0,353 1,348
3 -2,598 -2,488 1,481 -0,123 0,034 -1,57 -1,189 -2,793
4 0 7,294 4,242 0,04 6,063 1,861 5,793 1,591
5 -2,599 2,191 1,37 -0,233 0,444 -1,159 2,231 0,628
6 -5,197 -4,763 -0,814 0,181 -0,869 0,126 -2,817 -1,822
7 -2,597 5,261 1,838 0,233 2,25 0,645 2,835 1,23
8 0 3,83 3,956 -0,246 3,32 -0,882 4,431 0,229
9 -2,6 -1,645 1,312 -0,29 2,158 0,556 0,43 -1,172
10 -5,197 3,033 -0,703 0,292 -1,029 -0,034 -0,078 0,917
11 -2,596 0,293 1,947 0,341 1,309 -0,297 1,634 0,028
12 0 1,645 3,842 -0,36 4,519 0,317 3,466 -0,736
13 -2,601 5,79 1,2 -0,401 1,361 -0,24 2,258 0,657
14 -5,197 -3,179 -0,581 0,414 -0,887 0,108 -1,064 -0,069
15 -2,596 -0,126 2,051 0,445 1,696 0,09 1,165 -0,441
16 0 8,332 3,721 -0,481 3,958 -0,244 4,666 0,464
17 -2,601 -1,342 1,101 -0,5 1,74 0,139 1,484 -0,117
18 -5,197 -1,813 -0,46 0,535 -0,944 0,051 -1,256 -0,261
19 -2,595 5,463 2,144 0,537 1,528 -0,079 1,941 0,334
20 0 0,655 3,6 -0,602 4,181 -0,021 4,084 -0,118 
 
 
Hermitiano de ordem 3
-6
-4
-2
0
2
4
6
8
10
0 2 4 6 8 10 12 14 16 18 20
N (n. de passos)
x1
x1
x1 H3
R3
 
Fig. 8 - Aniquilamento no algoritmo de ordem 3 
 
 
 
 
Hermitiano de ordem 4
-6
-4
-2
0
2
4
6
8
10
0 2 4 6 8 10 12 14 16 18 20
N(n. passos)
x1
x1
x1 H4
R4
 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
24
 
Fig. 9 - Aniquilamento no algoritmo de ordem 4 
 
 
 
Hermitiano de ordem 5
-6
-4
-2
0
2
4
6
8
10
0 2 4 6 8 10 12 14 16 18 20
N (n. de passos)
x1
x1
x1 H5
R4
 
Fig. 10 - Aniquilamento no algoritmo de ordem 5 
 
 
11. CONCLUSÕES 
 
 Com o esgotamento dos algoritmos de ordens inferiores apontado pela 
bibliografia, um caminho para obter métodos de integração que contenham as principais 
propriedades espectrais é a formulação hermitiana livre na montagem de processos de 
ordens superiores, por apresentar congenitamente consistência e convergência, e 
facilmente permitir uma análise espectral que indique a presença ou não de estabilidade 
e amortecimento numérico, este último associado ao aniquilamento assintótico. 
 A geração de sistemas de equações de ordem mais elevada pode ser contornada 
com técnicas mais eficientes de solução, além de que maior precisão permite passos 
maiores. 
 Outra característica positiva dos métodos criados a partir dessa abordagem é a 
facilidade de introduzir parâmetros de ajuste. 
 A família apresentada é a priori incondicionalmente estável, com aniquilamento 
assintótico, cuja presença era um objetivo importante por razões acima descritas. 
Detalhou-se o caminho para obter-se a ordem de erro local que se desejar. 
 A regra de formação dos algoritmos está perfeitamente definida e é bastante clara. 
Aparentemente, dada uma ordem de precisão qualquer, é sempre possível gerar um 
método incondicionalmente estável e com aniquilamento assintótico presente. 
 O critério de tomar um operador bastante geral combinado com o decorrente da 
sua derivação em relação ao tempo e a subseqüente busca de requisitos para o 
atendimento de propriedades espectrais desejadas mostraram-se eficazes na 
identificação de algoritmos potencialmente viáveis. 
 Dentre os algoritmos deduzidos, observa-se que o primeiro coincide com o já 
obtido por Laier, em trabalho apresentado em 1995. 
 O algoritmo de segunda ordem resultante da presente abordagem utiliza também 
um operador empregado por Laier no trabalho citado, embora combinado 
diferentemente. Como resultado, obtém-se um método cujo desempenho espectral é 
idêntico ao de Laier, que, por sua vez, apresenta a mesma dissipação do método de 
Hulbert, não tão brusca quanto a do processo de primeira ordem da mesma família e do 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
25
método de Houbolt. O alongamento também se comporta como no método de Hulbert e 
de Laier. 
 Os algoritmos de ordens superiores têm desempenho espectral progressivamente 
melhor, como mostram os gráficos apresentados. 
 A dissipação surge para passos de tamanho crescente com a ordem do erro e de 
forma cada vez mais brusca, definindo cada vez mais nitidamente uma freqüência de 
corte. 
 A dispersão observada nos gráficos de alongamento versus tamanho do passo está 
praticamente ausente no início do gráfico, a partir da terceira ordem, porém aumenta 
abruptamente para ∆t/T > 0,47 aproximadamente. Esse valor não inviabiliza o uso do 
processo, sendo mesmo um limite superior ao de outros métodos existentes. 
 O exemplo numérico apresentado confirma e enfatiza as qualidades apontadas 
pela análise espectral. Para análise transiente, o ganho em rapidez nos métodos de 
ordens superiores é a principal vantagem, enquanto para respostas ao fim de intervalos 
de tempo grandes (long-term analysis) soma-se a isto o fato da resposta manter 
qualidade, tardando a apresentar a deterioração provocada pelo aniquilamento 
assintótico, característica que pode influir precocemente no resultado no caso das 
ordens inferiores. 
 Uma comparação mais precisa do tempo demandado no processamento versus a 
precisão obtida, embora relevante, seria muito particular para o caso de vibrações livres. 
Com a presença da excitação, os métodos de ordens superiores (e não apenas os aqui 
tratados) exigem operações numéricas envolvendo as funções que a representam, de 
modo que as conclusões também sejam relativas ao problema específico enfrentado. 
 Com a explicitação de algoritmos incondicionalmente estáveis e com o 
aniquilamento assintótico de ordens elevadas, abre-se a possibilidade de, introduzindo-
se parâmetros livres de ajuste às suas expressões, aprimorar-se uma ou mais de suas 
características espectrais, ou mesmo introduzir propriedades nas matrizes dos 
coeficientes que auxiliem a solução numérica. Essa é uma seqüência possívelpara o 
presente trabalho; outra que se vislumbra é a implementação dos algoritmos obtidos à 
solução de problemas específicos, visando confirmar suas potencialidades indicadas 
pela análise espectral. 
 
 
REFERÊNCIAS BIBLIOGRÁFICAS 
 
ARGYRIS, J.; DUNNE, P.C.; ANGELOPOULOS, L. Non-linear oscillations using 
the finite element technique. Comp. Meth. Appl. Mech. Eng., v.2, p. 203-250, 1973. 
ARGYRIS, J. H.; DOLTSINIS, J. St.; KNUDSON, W. C.; VAZ, L. E.; WILLAM, 
K. J. Numerical solution of transient nonlinear problems. Comp. Meth. Appl. Mech. 
Eng., v. 17/18, p. 341- 409, 1979. 
ARGYRIS, J.; MLEJNEK, H. P. Texts on computational mechanics: v.5 - Dynamics 
of structures. Amsterdam, North-Holland, 1991. 
ARGYRIS, J. H.; VAZ, L. E.; WILLAM, K. J. Higher order methods for transient 
diffusion analysis. Comp. Meth. Appl. Mech. Eng., v. 12, n. 2, p. 243-278, 1977. 
BATHE, K. J.; WILSON, E. L. Numerical methods in finite element analysis. 
Englewood Cliffs, Prentice-Hall, 1976. 
BAZZI, G.; ANDERHEGGEN, E. The ρ-family of algoithms for time-step integration 
with improved numerical dissipation. Earthqu. Eng. Struct. Dyn., v. 10, p. 537-
550, 1982. 
H.M. Bottura & J.E. Laier 
Cadernos de Engenharia de Estruturas, São Carlos, n.14, p. 1-27, 1999 
26
BRUSA, L.; NIGRO, L. A one-step method for direct integration os structural 
dynamic equations. Int. J. Numer. Meth. Eng., v. 15, p. 685-699, 1980. 
CHUNG, J.; LEE, J. M. A new family of explicit time integration methods for linear 
and non-linear structural dynamics. Int. J. Numer. Meth. Eng., v. 37, p. 3961-
3976, 1994. 
COLLATZ, L. The numerical treatment of differential equations. Berlin, Springer-
Verlag, 1960. 
GOUDREAU, G. L.; TAYLOR, R. L. Evaluation of numerical integration methods in 
elastodynamics. Comp. Meth. Appl. Mech. Eng., v. 2, p. 69-97, 1972. 
HILBER, H. M. ; HUGHES, T. J. R. Collocation, dissipation and overshoot for time 
integration schemes in structural dynamics. Earthqu. Eng. Struct. Dyn., v. 6, p. 99-
117, 1978. 
HOFF, C. ; PAHL, P. J. Development of an implicit method with numerical 
dissipation from a generalized single-step algorithm for structural dynamics. Comp. 
Meth. Appl. Mech. Eng., n. 67, p. 367-385, 1988. 
HOUBOLT, J. C. A recurrence matrix solution for the dynamic response of elastic 
aircraft. J. Aer. Sci., v. 17, p. 540-550, 1950. 
HUGHES, T. J. R. ; HULBERT, G. M. Space-time finite element methods for 
elastodynamics formulations and error estimates. Comp. Meth. Appl. Mech. Eng., v. 
66, p. 339-363, 1988. 
HULBERT, G. M. Time finite element methods for structural dynamics. Int. J. Numer. 
Meth. Eng., v. 33, p. 307-331, 1992. 
HULBERT, G. M. A unified set of single-step asymptotic annihillation algorithms for 
structural dynamics. Comp.Meth. Appl. Mech. Eng., v. 113, p. 1-9, 1994. 
HULBERT, G. M. Limitations on linear multistep methods for structural dynamics. 
Earthqu. Eng. Struct. Dyn., v. 20, p. 191-196, 1991. 
HULBERT, G. M.; CHUNG. J. The unimportance of the spurious root of time 
integration algorithms for structural dynamics. Commun. Numer. Meth. Eng., v. 
10, p. 591-597, 1994. 
HULBERT, G. M.; CHUNG, J. A family of single-step houbolt time integration 
algorithms for structural dynamics. Comp. Meth. Appl. Mech. Eng., v. 118, p. 1-
11, 1994. 
KATONA, M. G.; ZIENKIEWICZ, O. C. A unified set of single-step algorithms - 
Part 3: The Beta-M Method, a generalization of the newmark scheme. Int. J. 
Numer. Meth. Eng., v. 21, p. 1345-1359, 1985. 
KRIEG, R. D. Unconditional stability in numerical time integration methods. J. Appl. 
Mech., v. 40, p. 417-421, 1973. 
KUJAWSKI, J. ; GALLAGHER, R. H. A generalized least-squares family of 
algorithms for transient dynamic analysis. Earthqu. Eng. Struct. Dyn., v. 18, p. 
539-550, 1989. 
LAIER, J. E. Algoritmo hermitiano de integração passo a passo das equações da 
dinâmica das estruturas. In: JORNADAS SUDAMERICANAS DE INGENIERÍA 
ESTRUCTURAL, 26., Montevideo, Uruguay, 1993. Memorias. v. 1, p. 19-29. 
LAIER, J. E. Asymptotic annihilation hermitian algorithms for integration in time 
(submetido para publicação na revista Computer & Structures). 
MAKINSON, G. J. Stable high order implicit methods for the numerical solution of 
systems of differential equations. The Computer Journal, v. 11, n. 3, p. 305-310, 
1968. 
Uma família de algoritmos hermitianos para a integração direta das equações de dinâmica das estruturas 
Cadernos de Engenharia de Estruturas, São Carlos, n. 14, p. 1-27, 1999 
27
NEWMARK, N. M. A method of computation for structural dynamics. Proc. of ASCE, 
v. 85, n.EM3, p. 67-94, 1959. 
NØRSETT, S. P. One-step methods of hermite type for numerical integration of stiff 
systems. BIT : Nordisk Tidskriff for Inform. Behandung, n. 14, p. 63-77, 1974. 
PEZESHK, S.; CAMP, C. V. An explicit time integration technique for dynamic 
analyses. Int. J. Numer. Meth. Eng., v. 38, p. 2265-2281, 1995. 
WARBURTON, G. B. The Dynamical behaviour of structures. Oxford, Pergamon 
Press, 1976. 
WILSON, E. L. A computer program for the dynamic stress analysis of underground 
structures. Berkeley, Division of Struct. Eng. And Struct. Mech., Univ. of 
California, 1968. (SESM Report, N. 68-1). 
WOOD, W. L.; BOSSAK, M.; ZIENKIEWICZ, O. An alpha modification of 
Newmark's Method. Int. J. Numer. Meth. Eng., v. 15, p. 1562-1566, 1980. 
ZIENKIEWICZ, O. A new look at the newmark, houbolt and other time stepping 
formulas. a weighted residual approach. Earthqu. Eng. Struct. Dyn., v. 5, p. 413-
418, 1977. 
ZIENKIEWICZ, O. C.; WOOD, W. L.; TAYLOR, R. L. An alternative single-step 
algorithm for dynamic problems. Earthqu. Eng. Struct. Dyn., v. 3, p. 31-40, 1980.