MÉTODOS NUMÉRICOS E EQUAÇÕES DIFERENCIAIS EDO

Prévia do material em texto

MÉTODOS NUMÉRICOS EM 
EQUAÇÕES DIFERENCIAIS 
AULA 3 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Profª Marina Vargas Reis de Paula Gonçalves 
 
 
2 
EQUAÇÕES DE ORDEM SUPERIOR E SISTEMAS DE EQUAÇÕES 
DIFERENCIAIS 
Nas aulas anteriores, foram estudas equações diferenciais que envolviam 
somente a primeira derivada de uma função, ou seja, equações diferenciais de 
primeira ordem. No entanto, na modelagem de diversos problemas, surge a 
necessidade de utilizar derivadas de ordem superior. 
Dessa forma, passaremos a estudar a resolução de problemas regidos 
por equações diferenciais ordinárias de ordem superior, abrangendo modelos 
como sistemas massa-mola, deformações lineares de sólidos e circuitos 
elétricos RLC. 
Ao longo deste capítulo, serão apresentadas definições de EDOs de 
ordem superior e como transformá-las em uma análise de diversas EDOs de 
primeira ordem, acopladas por meio de um sistema de equações. 
Com isso, serão revisitados alguns métodos para solução numérica 
apresentados anteriormente, com o enfoque de solução de sistemas de 
equações diferenciais ordinárias de primeira ordem. Nesse contexto, serão 
apresentados aspectos que envolvem a estabilidade de solução numérica no 
tempo, bem como características inerentes de certas EDOs, chamadas rígidas, 
que devem ser observadas para a sua correta modelagem numérica. 
Os assuntos abordados no escopo deste material podem ser encontrados 
nos materiais de Burden e Faires (2008), Franco (2006), Butcher e Goodwin 
(2008) e Chapra (2012). Esses materiais são indicados para um maior 
aprofundamento do conteúdo. 
TEMA 1 – EDO DE ORDEM SUPERIOR E O SISTEMA DE EQUAÇÕES 
DIFERENCIAIS ASSOCIADO 
Apesar de a extensão dos conceitos ser poderosa no sentido de 
aplicabilidade, a estratégia numérica de solução herda vários aspectos já 
discutidos para o caso de primeira ordem. Isso se deve ao fato de que iremos 
procurar traduzir a equação diferencial de ordem superior em um sistema de 
equações diferenciais de primeira ordem por meio de substituições de variáveis 
simples. 
 
 
 
3 
1.1 Substituição de variáveis para gerar o sistema de EDOs 
Inicialmente, estenderemos a notação do problema de valor inicial 
apresentada anteriormente como: 
𝒅𝒖
𝒅𝒕
= 𝒇(𝒕, 𝒖(𝒕)) 
Desse modo, contemos derivadas de ordem 𝒏. Assim, escrevemos: 
𝒅𝒏𝒖
𝒅𝒕𝒏
= 𝒇 (𝒕, 𝒖,
𝒅𝒖
𝒅𝒕
,
𝒅𝟐𝒖
𝒅𝒕𝟐
, … ,
𝒅𝒏−𝟏𝒖
𝒅𝒕𝒏−𝟏
) 
Define-se assim uma equação diferencial de ordem superior de maneira 
análoga. 
Naturalmente, é necessário estipular condições iniciais para cada uma 
das derivadas até a ordem 𝒏 − 𝟏. Ou seja, 𝒖(𝒕𝟎) = 𝜷𝟎
𝟎, 
𝒅𝒖
𝒅𝒕
(𝒕𝟎) = 𝜷𝟎
𝟏, 
𝒅𝟐𝒖
𝒅𝒕𝟐
(𝒕𝟎) =
𝜷𝟎
𝟐, 
𝒅𝟑𝒖
𝒅𝒕𝟑
(𝒕𝟎) = 𝜷𝟎
𝟑, …, 
𝒅𝒏−𝟏𝒖
𝒅𝒕𝒏−𝟏
(𝒕𝟎) = 𝜷𝟎
𝒏−𝟏. 
Dessa forma, podemos transformar a equação diferencial em um sistema 
de equações diferenciais associando cada derivada de ordem 𝒋, com 𝒋 de 𝟎 a 
𝒏 − 𝟏, a uma variável 𝒗𝒋 =
𝒅𝒋𝒖
𝒅𝒕𝒋
 e, assim, a EDO fica: 
𝒅𝒏𝒖
𝒅𝒕𝒏
= 𝒇(𝒕, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏) 
Com as respectivas condições iniciais. É importante notar que, ao 
escolhermos 𝒗𝒋 =
𝒅𝒋𝒖
𝒅𝒕𝒋
, construímos um sistema de equações diferenciais 
lineares. 
TEMA 2 – SOLUÇÃO USANDO MÉTODOS DE EULER 
Os métodos apresentados de passo simples para solução numérica de 
um problema de valor inicial (PVI) podem ser facilmente estendidos para a 
solução de sistemas de equações diferenciais sujeita à imposição de um valor 
inicial. 
Esse procedimento consiste em simplesmente dar um passo no tempo 
utilizando o método escolhido para cada equação do sistema. 
 
 
 
 
 
4 
2.1 Método explícito de Euler 
Para o caso do método explícito de Euler, escrevemos o processo iterativo 
para cada uma das 𝒏 equações diferenciais provenientes da EDO de ordem 𝒏. 
Ou seja, para cada uma das equações, temos: 
𝒅𝒗𝒋
𝒅𝒕
= 𝒇(𝒕, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏) 
Portanto: 
𝒗𝒋(𝒕
𝒊+𝟏) = 𝒗𝒋(𝒕
𝒊) + 𝒉𝒇(𝒕, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏) 
2.2 Método de Euler modificado 
O método de Euler apresentado pode ser repensado para ter uma 
aproximação melhor, ao se utilizar uma função 𝒇 constante no intervalo de 
integração temporal, mas sim como uma média entre o valor inicial e final. 
Ou seja: 
�̅� =
𝒇𝒊+𝟏(𝒕, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏) + 𝒇𝒊(𝒕, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏)
𝟐
 
Assim, 𝒇𝒊 denota a função 𝒇 calculada no passo de tempo 𝒊. 
Os métodos para resolver as equações diferenciais de primeira ordem são 
generalização dos métodos para equações de primeira ordem. Para aproximar 
a solução de um problema de valor inicial de segunda ordem: 
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥, 𝑦, 𝑦´), 𝑦(𝑥0) = 𝑦0, 𝑦´(𝑥0) = 𝑦´0 
Reduzimos a equação diferencial de segunda ordem a um sistema de 
primeira ordem. Fazendo 𝑦´ = 𝑢, a equação se torna: 
{
𝑦´ = 𝑢
𝑢´ = 𝑓(𝑥, 𝑦, 𝑢)
 
Aplicamos em seguida um método particular para cada equação do 
sistema resultante. Por exemplo, o método de Euler seria: 
{
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑢𝑛
𝑢𝑛+1 = 𝑢𝑛 + ℎ𝑓(𝑥𝑛, 𝑦𝑛, 𝑢𝑛)
. 
 
 
 
 
5 
TEMA 3 – SOLUÇÃO USANDO O MÉTODO DE RUNGE-KUTTA DE QUARTA 
ORDEM 
De forma análoga ao que foi desenvolvido para o método de Euler no 
sistema de equações diferenciais decorrente da EDO da ordem superior, 
podemos estender a aplicação do método de Runge-Kutta de quarta ordem para 
esses casos. 
3.1 As constantes de taxa de variação instantânea 
No método de Runge-Kutta de quarta ordem, a aproximação se dá de 
forma iterativa, com base em uma inclinação da predição calculada com base 
em quatro taxas de variação auxiliares (𝒌𝟏, 𝒌𝟐, 𝒌𝟑, 𝒌𝟒). 
Primeiramente, calcula-se 𝒌𝟏em relação ao ponto inicial do intervalo de 
integração no tempo (𝒊). Com base em 𝒌𝟏, é calculada a inclinação no ponto 
médio do intervalo (𝒕𝒊 + 𝒉/𝟐), 𝒌𝟐. Com essas duas constantes, retorna-se ao 
ponto inicial do intervalo para estimar novamente a taxa de variação no ponto 
médio, 𝒌𝟑. Por fim, com 𝒌𝟏, 𝒌𝟐 e 𝒌𝟑, calcula-se 𝒌𝟒, referente ao ponto final do 
intervalo de integração no tempo (𝒕𝒊+𝟏). 
Para aproximar a solução de um problema de valor inicial de segunda 
ordem: 
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥, 𝑦, 𝑦´), 𝑦(𝑥0) = 𝑦0, 𝑦´(𝑥0) = 𝑦´0 
Da mesma forma que, para o método de Euler, reduzimos a equação 
diferencial de segunda ordem a um sistema de primeira ordem. Fazendo 𝑦´ = 𝑢, 
a equação se torna: 
{
𝑦´ = 𝑢
𝑢´ = 𝑓(𝑥, 𝑦, 𝑢)
 
Aplicamos em seguida um método particular para cada equação do 
sistema resultante, o Runge-Kutta de quarta ordem. 
Dessa forma, as constantes 𝒌 são calculadas para o caso de múltiplas 
EDOs, como: 
𝒌𝟏 = 𝒇 (𝒕
𝒊, 𝒖(𝒕𝒊)) = 𝒇(𝒕𝒊, 𝒗𝟎, 𝒗𝟏, 𝒗𝟐, … , 𝒗𝒏−𝟏) 
𝒌𝟐 = 𝒇 (𝒕
𝒊 +
𝒉
𝟐
, 𝒗𝟎(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟏, 𝒗𝟏(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟏, 𝒗𝟐(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟏, … , 𝒗𝒏−𝟏(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟏) 
 
 
6 
𝒌𝟑 = 𝒇 (𝒕
𝒊 +
𝒉
𝟐
, 𝒗𝟎(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟐, 𝒗𝟏(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟐, 𝒗𝟐(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟐, … , 𝒗𝒏−𝟏(𝒕
𝒊) +
𝒉
𝟐
𝒌𝟐) 
𝒌𝟒 = 𝒇(𝒕
𝒊+𝟏, 𝒗𝟎(𝒕
𝒊) + 𝒌𝟑, 𝒗𝟏(𝒕
𝒊) + 𝒌𝟑, 𝒗𝟐(𝒕
𝒊) + 𝒌𝟑, … , 𝒗𝒏−𝟏(𝒕
𝒊) + 𝒌𝟑) 
3.2 Procedimento iterativo 
De posse das constantes 𝒌, o procedimento iterativo é facilmente descrito 
por meio de uma taxa média ponderada de variação no intervalo por meio da 
seguinte expressão: 
𝒗𝒋(𝒕
𝒊+𝟏) = 𝒗𝒋(𝒕
𝒊) +
𝒉
𝟔
[𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒] 
TEMA 4 – ESTABILIDADE 
A escolha de determinado método numérico permeia diversos aspectos, 
como grau de acurácia desejado, custo computacional e tipo de aplicado. Além 
desses aspectos, há casos em que determinados métodos têm garantias 
teóricas de que funcionarão como desejado, enquanto outros não são 
adequados. 
Nesse sentido, é preciso avaliar os métodos de escolha conforme sua 
estabilidade em determinada aplicação. 
No contexto de estabilidade numérica, há dois conceitos-chave: 
convergência e consistência. Esses conceitos serão definidoscom precisão a 
seguir. Por fim, é apresentado um resultado que define de forma objetiva a 
estabilidade de um método numérico no contexto da solução numérica de 
problemas de valor inicial. 
4.1 Convergência e consistência 
Primeiramente, dizemos que um método é convergente se o refino 
arbitrário da malha implica em aproximação da resposta aproximada em relação 
à resposta verdadeira do problema. Ou seja, um método é convergente se: 
𝐥𝐢𝐦
𝒉→𝟎
𝐦𝐚𝐱
𝟎≤𝒊≤𝑵−𝟏
|�̅�𝒊 − 𝒖𝒊| = 𝟎 
Em que �̅�𝒊 aproxima 𝒖𝒊 no 𝒊-ésimo passo em uma malha temporal de 
abertura 𝒉. 
Por outro lado, por se tratar de uma aproximação numérica, é necessário 
se preocupar em como os erros de truncamento impactam a solução. Assim, um 
 
 
7 
método é dito consistente se o erro de truncamento 𝝉(𝒉) associado a uma malha 
de abertura 𝒉 diminui arbitrariamente à medida que a malha é refinada. Ou seja: 
𝐥𝐢𝐦
𝒉→𝟎
𝐦𝐚𝐱
𝟎≤𝒊≤𝑵−𝟏
|𝝉𝒊(𝒉) | = 𝟎 
4.2 Estabilidade de um método de PVI 
Consideremos um método numérico para o problema de valor inicial de 
primeira ordem, escrito de forma que a solução aproximada �̅� seja obtida como: 
�̅�𝒊+𝟏 = �̅�𝒊 + 𝒉 ∗ 𝒈(𝒕𝒊, 𝒖𝒊, 𝒉) 
Suponha que exista um tamanho máximo de malha 𝒉𝟎 > 𝟎, tal que 
𝒈(𝒕𝒊, 𝒖𝒊, 𝒉) seja Lipschitz contínua em relação à variável 𝒖 com constante de 
Lipschitz de valor 𝒌 para o domínio: 
𝑫 = {𝒈(𝒕, 𝒖, 𝒉)/ 𝐭𝒊𝒏𝒊𝒄𝒊𝒂𝒍 ≤ 𝐭 ≤ 𝐭𝒇𝒊𝒏𝒂𝒍, 𝒖 ∈ 𝕽, 𝟎 ≤ 𝒉 ≤ 𝒉𝟎} 
Então podemos afirmar que: 
 O método é estável; 
 O método é consistente se, e somente se, é convergente; 
 Se o erro de truncamento em um nó da malha |𝝉𝒊(𝒉)| é limitado 
superiormente pelo erro de truncamento global 𝝉(𝒉) dentro de um 
intervalo de malha 𝟎 ≤ 𝒉 ≤ 𝒉𝟎, então: 
|�̅�𝒊 − 𝒖𝒊| ≤
𝝉(𝒉)
𝒌
𝒆𝒌(𝒕−𝒕𝒊𝒏𝒊𝒄𝒊𝒂𝒍) 
TEMA 5 – EQUAÇÕES DIFERENCIAIS RÍGIDAS 
Podemos entender equações diferenciais rígidas como aquelas em que 
termos de grande variação ao longo do tempo coexistem com termos de pouca 
variação. Esse tipo de equação diferencial se caracteriza intrinsecamente pela 
natureza de variação abrupta na resposta do problema e/ou em suas derivadas. 
Em geral, problemas rígidos apresentam solução analítica contendo um 
termo da forma 𝒆−𝒔𝒕 com 𝒔 ≫ 𝟏. Dessa forma, esse termo tende a sumir 
rapidamente à medida que se avança no tempo. No entanto, seu efeito na 𝒏-
ésima derivada é da forma 𝒔𝒏𝒆−𝒔𝒕, a qual tende a persistir, dada a magnitude de 
𝒔 e, portanto, há uma dificuldade em manter limitada a respectiva derivada, 
impactando na estabilidade da solução numérica. 
 
 
 
8 
5.1 Condições de estabilidade no tempo 
Analisemos inicialmente uma EDO de primeira ordem, como as que 
aparecem em cada uma das equações diferenciais do sistema gerado por uma 
equação diferencial de ordem superior. Por exemplo, temos uma EDO da forma: 
𝒅𝒖
𝒅𝒕
= −𝒌𝒖 
Essa forma está sujeita à condição inicial 𝒖(𝒕 = 𝟎) = 𝒖𝟎. Essa EDO tem 
como solução analítica a expressão: 
𝒖 = 𝒖𝟎𝒆
−𝒌𝒕 
Se formos resolver essa EDO numericamente pelo método de Euler, por 
exemplo, teremos o seguinte procedimento iterativo: 
𝒖𝒊+𝟏 = 𝒖𝒊 + 𝒉
𝒅𝒖
𝒅𝒕
 
Ou seja, utilizando o conhecimento da EDO em questão, temos: 
𝒖𝒊+𝟏 = 𝒖𝒊 − 𝒌𝒖𝒊𝒉 
𝒖𝒊+𝟏 = 𝒖𝒊(𝟏 − 𝒌𝒉) 
Considerando que desejamos que a variável temporal em questão seja 
limitada no tempo, devemos impor que o fator multiplicativo (𝟏 − 𝒌𝒉) contenha a 
evolução temporal. Dessa forma, devemos ter |𝟏 − 𝒌𝒉| < 𝟏. Portanto, para 
termos uma resposta estável no tempo, devemos escolher 𝒉 limitado, de forma 
que 𝒉 <
𝟐
𝒌
. Nota-se nesse caso que o método explícito de Euler, como formulado, 
é condicionalmente estável. 
NA PRÁTICA 
Considere o PVI associado ao modelo oscilação de um pêndulo simples 
de massa 𝒎, comprimento de corda 𝒍, de posição angular em relação à vertical 
𝜽, posição inicial 𝜽(𝟎) = 𝜽𝟎, velocidade angular inicial de 
𝒅𝜽
𝒅𝒕
(𝟎) = �̇�𝟎 e 
aceleração da gravidade 𝒈. 
Com essa configuração física, o PVI pode ser escrito com base no 
equilíbrio energético: 
𝒎𝒍
𝒅𝟐𝜽
𝒅𝒕𝟐
= 𝒎𝒈𝜽 
Deve-se considerar que a 𝜽 varia em valores baixos. Dessa forma, 
podemos escrever: 
 
 
9 
𝒅𝟐𝜽
𝒅𝒕𝟐
−
𝒈
𝒍
𝜽 = 𝟎 
Adote 𝜽(𝟎) = 𝟎, 
𝒅𝜽
𝒅𝒕
(𝟎) = 𝟎, 𝟓 𝒓𝒂𝒅/𝒔, 𝒍 = 𝟎, 𝟒 𝒎, 𝒈 = 𝟗, 𝟖𝟏
𝒎
𝒔𝟐
 e um intervalo 
de tempo de análise de acordo com a necessidade para observar com 
comportamento da resposta. Com base nessa contextualização: 
 Faça uma análise da estabilidade para o método de Euler explícito e 
determine um valor crítico de 𝒉, se houver; 
 Utilize o método de Euler para calcular a resposta ao longo do tempo para 
𝜽(𝒕) e 
𝒅𝜽
𝒅𝒕
(𝒕) para dois valores distintos de 𝒉, um pequeno e um grande, e 
discuta a diferença nas respostas; 
 Utilize o método de Euler modificado para calcular a resposta ao longo do 
tempo para 𝜽(𝒕) e 
𝒅𝜽
𝒅𝒕
(𝒕) para dois valores distintos de 𝒉, um pequeno e 
um grande, e discuta a diferença nas respostas; 
 Utilize o método de Runge-Kutta de quarta ordem para calcular a resposta 
ao longo do tempo para 𝜽(𝒕) e 
𝒅𝜽
𝒅𝒕
(𝒕) para dois valores distintos de 𝒉, um 
pequeno e um grande, e discuta a diferença nas respostas; 
 Por fim, analise os resultados e discuta qual dos métodos você escolheria 
para a modelagem numérica do problema, considerando respostas para 
𝜽(𝒕) e 
𝒅𝜽
𝒅𝒕
(𝒕). 
Uma observação interessante diz respeito a uma simplificação desse 
modelo que possibilita a construção de uma EDO linear. Essa consideração 
consiste em admitir que o ângulo possui uma magnitude suficiente pequena. No 
caso geral, devemos considerar a altura real que afeta diretamente a energia 
potencial gravitacional. Assim, a EDO admite uma forma não linear, como: 
𝒅𝟐𝜽
𝒅𝒕𝟐
−
𝒈
𝒍
𝒔𝒆𝒏(𝜽) = 𝟎 
Esse tipo de problema será objeto de estudo futuramente. Por enquanto, 
é suficiente considerarmos o caso simplificado e linear. 
 
 
 
10 
FINALIZANDO 
Foram reapresentados métodos para a solução numérica das seções 
anteriores com o enfoque de solução de sistemas de equações diferenciais 
ordinárias de primeira ordem. Nesse contexto, foram discutidos aspectos que 
envolvem a estabilidade de solução numérica no tempo e as características 
inerentes de certas EDOs chamadas rígidas, que devem ser observadas para a 
sua correta modelagem numérica. 
Estudou-se como a escolha do passo 𝒉 pode interferir diretamente na 
estabilidade da solução, principalmente para equações diferenciais rígidas. 
 
 
 
11 
REFERÊNCIAS 
BURDEN, R. L.; FAIRES, J. D. Análise numérica. São Paulo: Cengage 
Learning, 2008. 
CHAPRA, S. C. Applied numerical methods. Columbus: McGraw-Hill, 2012. 
BUTCHER, J. C.; GOODWIN, N. Numerical methods for ordinary differential 
equations. New York: Wiley, 2008. 
FRANCO, N. B. Cálculo numérico. São Paulo: Pearson, 2006.