Uma equação diferencial ordinária —ou de forma abreviada, EDO— de ordem é uma relação entre uma função e as suas derivadas , , …, . Por exemplo, uma possível equação diferencial ordinária de terceira ordem é:
Neste caso, é possível também escrever a equação mostrando explicitamente a expressão que define a derivada de terceira ordem:
É conveniente usar uma notação que permite ver a equação de forma mais geral: a função será denotada por , será a primeira derivada e a segunda derivada . A equação diferencial é então:
que define a derivada de em relação a ( , , , ). Para resolver este problema é necessário resolver simultaneamente as equações que definem as derivadas de e . Ou seja, a equação original 3.1, de terceira ordem, foi transformada no sistema de 3 equações:
Esta forma de escrever a equação diferencial tem também a vantagem de poder ser escrita de forma mais compacta; define-se o vetor,
e o sistema de 3 equações diferenciais é equivalente a uma única equação diferencial vetorial:
onde a derivada do vetor é outro vetor obtido derivando cada uma das suas componentes,
e as três componentes da função vetorial são, neste caso,
Como tal, para resolver equações diferencias ordinárias como, por exemplo, a equação 3.1, basta saber resolver equações de primeira ordem com a forma geral
e generalizar o método ao caso em que e são vetores com várias componentes.
A forma geral das EDO de primeira ordem é
A função chama-se variável dependente e é a variável independente. Uma solução da EDO, num intervalo [ , ], é qualquer função de que verifica a equação.
Existem em geral muitas soluções, por exemplo, a figura 3.1 mostra 7 possíveis soluções da equação . Em cada ponto de cada uma das curvas, o declive tem o mesmo valor obtido substituindo as coordenadas do ponto na expressão . As diferentes soluções não se cruzam nunca, porque em cada ponto existe apenas um valor e cada ponto pertence unicamente a uma das soluções da equação diferencial.
Se for dado um valor inicial para a função no ponto inicial , existe uma única solução , que é a curva tal que e com declive igual a em qualquer valor de .
Os métodos de resolução numérica consistem em calcular o valor da variável dependente numa sequência discreta de pontos { , , , …, }, usando alguma aproximação. O resultado obtido aplicando um determinado método será o conjunto de pontos {( , ), ( , ), …, ( , )}, que aproximam o gráfico da função . Nos métodos que usam incrementos constantes, o intervalo [ , ] divide-se em subintervalos de comprimento , de forma que cada valor na sequência { } é igual ao valor anterior mais :
Por exemplo, a figura 3.2 mostra a solução da equação diferencial com e condição inicial = 1, no intervalo 0≤ ≤ 4. Os cinco pontos sobre a curva são a aproximação da função usando apenas 5 pontos com abcissas 0, 1, 2, 3 e 4.
Cada ângulo na figura é a inclinação segmento entre os pontos ( , ) e ( , ). A tangente do ângulo (declive médio) é igual ao valor médio de no intervalo [ , ]. A partir das coordenadas ( , ) de um dos pontos, o declive médio permite calcular as coordenadas do ponto seguinte:
Usando estas relações de recorrência de forma iterativa, consegue-se obter as coordenadas de todos os pontos ( , ), começando com os valores iniciais dados ( , ).
Se os declives médios pudessem ser calculados de forma exata, os resultados obtidos seriam exatos. No entanto, para calcular o valor médio de num intervalo seria necessário conhecer a função , mas quando se usam métodos numéricos é porque não existem métodos para encontrar a expressão analítica dessa função. Os diferentes métodos numéricos que existem para resolver equações diferenciais correspondem a diferentes esquemas para estimar o valor médio aproximado da função em cada intervalo.
No método de Euler admite-se que ≈ . Ou seja, o valor médio de em cada intervalo é aproximado pelo valor de no ponto inicial do intervalo. Realizando os cálculos no Maxima para o caso considerado acima, em que e com condição inicial = 1, pode usar-se uma lista para armazenar as coordenadas dos pontos; inicia-se a lista com o ponto inicial e define-se a função dada:
Os quatro pontos seguintes na sequência são acrescentados à lista , usando um ciclo e as relações de recorrência 3.12 com = 1:
O gráfico na figura 3.3 mostra o resultado obtido (segmentos de reta e pontos), comparando-o com a solução exata (curva contínua). Observa-se que a solução obtida com o método de Euler não é uma aproximação muito boa. No entanto, reduzindo o valor dos incrementos em de = 1 para um valor menor, deverá ser possível obter uma aproximação melhor. É conveniente reduzir gradualmente o valor de e comparar as soluções com as obtidas com valores maiores de ; se a solução obtida não variar significativamente quando é reduzido, essa solução deverá ser uma boa aproximação da solução verdadeira.
A carga armazenada num condensador ligado a uma pilha e a uma resistência é uma função que depende do tempo e verifica a seguinte equação diferencial
No instante inicial, = 0, a carga no condensador é nula. Usando o método de Euler, encontre uma sequência que aproxime no intervalo 0 ≤ ≤ 6.
Resolução. Neste caso, a função que define a derivada não depende da variável independente . Os valores iniciais e a função que define a derivada são:
Começando com incrementos de tempo =0.1, são necessárias 60 iterações para terminar em = 6.
O último ponto na lista é:
Convém repetir o processo, com um valor menor de , por exemplo 0.01 e comparar o novo resultado com o resultado anterior:
A pesar de que o resultado final coincide em 4 algarismos significativos nos dois casos, é interessante comparar as duas sequências completas. As duas soluções obtidas com = 0.1 e = 0.01 estão armazenadas nas listas e e podem ser representadas num gráfico. A figura 3.4 mostra que, a pesar de os resultados obtidos para = 6 serem quase idênticos nos dois casos, há uma discrepância entre as duas soluções obtidas, mais visível no intervalo 0.5 < < 3.5. No entanto, como essa discrepância não é muito grande, espera-se que a solução com = 0.01 já esteja bastante próxima da solução verdadeira. O gráfico foi feito com o seguinte comando
A maior discrepância entre a sequência de pontos que aproximam a solução exata na figura 3.2 e a sequência obtida usando o método de Euler, apresentada na figura 3.3, é no segundo ponto, ( , ). O valor médio do declive no primeiro intervalo, que era aproximadamente igual a zero, no método de Euler foi substituído por 2, que é o declive no ponto inicial ( , ). A figura 3.2 mostra que, a pesar de que os dois primeiros valores da função, e , são iguais, o declive muda de um valor positivo em = 0 para um valor negativo em = 1. No método de Euler usou-se o declive em como valor médio o valor do declive nesse intervalo. Se o valor médio do declive fosse aproximado pela média entre os dois declives em e , o resultado seria muito melhor.
O método de Euler melhorado consiste em admitir que o declive médio é igual à média entre os declives no ponto inicial ( , ) e no ponto final ( , ); ou seja,
O problema com esta equação é que para poder-se calcular seria necessário saber o valor de , mas esse valor é desconhecido enquanto a solução da equação não for conhecida. É possível fazer uma estimativa inicial do valor que poderá ter e a seguir melhorar essa estimativa. Uma primeira estimativa muito rudimentar consiste em dizer que é igual a .
Com essa abordagem, a resolução da mesma equação considerada na secção anterior, com valor inicial = 1 e com 5 pontos no intervalo 0 ≤ ≤ 4, pode ser feita assim: começa-se por iniciar a lista de pontos e definir a função que calcula o declive
Foi necessário definir novamente a função , já que a definição usada na secção anterior foi logo substituída pela função do exemplo 3.1. Basta repetir as 4 iterações mas usando agora a expressão 3.13 nas relações de recorrência 3.12
A figura 3.5 mostra a sequência de pontos obtida, juntamente com a solução exata (curva contínua na figura). O resultado é muito melhor que o resultado obtido com o método de Euler.
O resultado obtido com (%i15) pode ser melhorado se a estimativa inicial = for melhorada com o resultado obtido a partir dessa primeira estimativa. Se designarmos à primeira estimativa do valor de , então = , como foi usado em (%i15). A segunda estimativa é: . Esse valor estará mais próximo do valor real de do que a estimativa inicial ; assim sendo, espera-se que a sequência se aproxime do valor real de . Calculam-se alguns termos dessa sequência, até ser menor que uma precisão desejada. Nesse momento usa-se = como valor inicial para o seguinte intervalo e o processo repete-se iterativamente em todo o intervalo de .
No caso do mesmo problema resolvido em (%i15), com = 1 e usando precisão de 0.0001 para , o método pode ser implementado assim:
Definiu-se novamente a função , acrescentando o comando float, para evitar obter fracções com numeradores e denominadores muito grandes; nos exemplos anteriores foram feitos menos cálculos e, por isso, não produziram frações com números grandes. O resultado (ver figura 3.6) mostra que a solução obtida está mais próxima da verdadeira solução do que as aproximações apresentadas nas figuras 3.3 e 3.5.
A pesar de se ter usado uma precisão de 0.0001, não se pode garantir que a solução obtida esteja dentro dessa precisão em relação á solução real, já que a equação 3.13 é apenas uma aproximação e não o valor médio real da derivada. Para melhorar essa aproximação, em vez de se obter a média do declive em dois pontos será necessário fazer uma média em vários pontos do intervalo em questão, com diferentes pesos que podem ser definidos de forma a minimizar o erro. Existem muitos métodos diferentes que usam esse esquema; o mais popular é o que será descrito na seguinte secção.
Neste método o valor médio do declive obtém-se a partir da média dos declives em 4 pontos diferentes, com pesos diferentes (figura 3.7).
Começa-se por calcular o declive no ponto inicial do intervalo, tal como no método de Euler:
a seguir, realiza-se um deslocamento na direção desse declive, avançando-se uma distância no eixo das abcissas, até um ponto 2 (ver figura 3.7). Nesse ponto 2, calcula-se um segundo valor do declive
Esse novo valor do declive é usado novamente, para realizar outro deslocamento a partir do ponto inicial, avançando na direção do eixo das abcissas, até um outro ponto 3, onde é calculado um terceiro valor do declive
seguindo a direção desse declive , realiza-se um terceiro deslocamento, a partir do ponto inicial, desta vez avançando uma distância no sentido do eixo das abcissas, para chegar até um ponto 4, onde se calcula um quarto valor do declive
Pode mostrar-se que para minimizar o erro cometido, o valor médio do declive deve ser aproximado pela seguinte combinação linear dos quatro declives calculados:
no exemplo da figura 3.7, esse valor médio da derivada desloca o ponto inicial até o ponto 4, que está bastante próximo da solução exata da equação.
Para aplicar este método à mesma equação considerada nas secções anteriores, com valor inicial = 1 e com 5 pontos no intervalo 0 ≤ ≤ 4, inicia-se a lista de pontos (a função que define o declive em cada ponto é a função que já está definida)
Neste caso foi dado explicitamente o valor do comprimento dos intervalos, para que as expressões seguintes sejam mais claras.
A seguir realizam-se as 4 iterações necessárias para chegar até o ponto = 4. Em cada iteração calculam-se os quatro declives , , e e a média com os pesos dados na equação 3.18.
A figura 3.8 mostra a sequência de pontos obtida, juntamente com a solução exata (curva contínua na figura). O resultado é bastante bom, a pesar do valor tão elevado que foi usado para os incrementos da variável . No caso da equação diferencial resolvida neste exemplo, existem métodos analíticos que permitem encontrar a expressão para a curva contínua que foi apresentada na figura; nos casos em que não é possível encontrar a solução analítica, o método de Runge-Kutta de quarta ordem é uma boa opção para encontrar uma boa aproximação numérica.
Já existe uma função rk predefinida no Maxima que executa o método de Runge-Kutta de quarta ordem. No exemplo acima, bastava indicar a expressão de , o nome que identifica a variável dependente, o seu valor inicial e o intervalo em que se quer obter a solução, incluindo o valor dos incrementos . Assim sendo, a mesma lista obtida com os comandos (%i17) até (%i21) podia ser obtida com um único comando:
Os métodos descritos nas secções anteriores podem ser generalizados facilmente a um sistema de várias equações, quando sejam conhecidas condições iniciais para todas as variáveis. Tal como se explica no inicio do capítulo, basta admitir que é um vetor e que é uma função vetorial.
Usando o método de Runge-Kutta de quarta ordem, encontre a solução do seguinte sistema de equações diferenciais:
no intervalo, 0 ≤ ≤ 4, com condições iniciais = -1.25 e = 0.75, em = 0.
Resolução. A variável independente é e a variável dependente é o vetor , que no Maxima representa-se por uma lista com duas partes, [ , ]. A expressão que define as derivadas das variáveis dependentes também será uma lista com duas partes, [ , ], mas é mais conveniente (para poder generalizar ao caso de variáveis) representá-la em função das componentes da lista :
foi dada também a variável independente como argumento da função, a pesar de que neste caso não depende de , para que os comandos usados a seguir possam ser reutilizados em casos mais gerais.
Com as condições iniciais = 0 e = (-1.25, 0.75) cria-se a lista onde serão guardados os resultados das iterações
Observe-se que a lista inicial [0,-1.25,0.75] e todas as outras listas que serão acrescentadas, incluem a variável independente e o vetor dependente . Para extrair apenas a variável independente usa-se o comando first, que extrai o primeiro elemento de uma lista e para extrair o vector usa-se o comando rest, que produz uma nova lista excluindo o primeiro elemento da lista original; ou seja, rest([0, -1.25, 0.75]) produz [-1.25, 0.75].
Usando incrementos = 0.02 para a variável independente , será necessário realizar 200 iterações para terminar em = 4.
O comando cons é semelhante a endcons, mas insere um elemento no início de uma lista, em vez de no fim. Neste caso, cons foi usado para inserir o valor de na lista [ , ], produzindo a lista [ , , ], que foi logo acrescentada no fim da lista dos resultados.
O último ponto na solução obtida é:
O mesmo resultado podia ter sido obtido com o comando rk
Para imprimir os resultados é conveniente usar printf. Esse comando aceita um formato do Lisp que diz como devem ser imprimidos os dados. Neste exemplo pode usar se o formato "∼{∼{∼f∼^, }∼%∼}" para imprimir os 3 elementos , e de cada iteração numa linha separada, separados por uma vírgula e um espaço. Os dois símbolos ∼{ e ∼} usam-se para delimitar um formato que deve ser aplicado recursivamente a todos os elementos de uma lista, neste caso a lista . Dentro desses símbolos, o formato ∼{∼f∼^, ∼}∼% aplicará recursivamente o formato ∼f ^, a cada elemento das sublistas de , neste caso a , e , seguido de um de fim de linha, representado pelo símbolo ∼%. O símbolo ∼f é usado para escrever números de vírgula flutuante e o símbolo ∼^ indica que o que vem a seguir só será escrito se não se tiver chegado ao último elemento da lista; ou seja, será escrita uma vírgula e um espaço após e , mas não após .
Para imprimir unicamente os 9 resultados para = 0, 0.5,… ,4, extraem-se unicamente esses elementos da lista , usando o comando makelist:
Finalmente, pode ser conveniente gravar esses resultados num ficheiro. Por exemplo, os programas de folha de cálculo conseguem ler dados como estão apresentados acima, dentro de ficheiros do tipo CSV (Comma Separated Values). Para gravar o que é apresentado no ecrã num ficheiro, usa-se o comando with_stdout
Os resultados ficam gravados no ficheiro resultados.csv