RUNGE E KUTTA DE 4ª ORDEM E MÉTODO DE ADAMS-BASHFORTH-MOULTON DE 4ª ORDEM

Prévia do material em texto

Universidade Federal do Maranhão - UFMA 
Bacharelado Interdisciplinar em Ciência e Tecnologia 
 
Camila Ferreira dos Santos 
Mariana Andressa Luna Pinheiro 
 
Professor Dr Kênio 
Cálculo Numérico 
Terceiro Trabalho Prático 
 
 
São Luís - MA 
2016 
 
1 
 
SUMÁRIO 
1. INTRODUÇÃO..............................................................................................................2 
2. FUNDAMENTAÇÃO TEÓRICA .................................................................................2 
3. PROBLEMA DE VALOR INICIAL..............................................................................4 
3.1 RUNGE E KUTTA DE 4ª ORDEM ..............................................................................5 
3.1.1 Resolução e implementação .............................................................................5 
3.2 MÉTODO DE ADAMS-BASHFORTH-MOULTON DE 4ª ORDEM ..........................7 
 3.2.1. Resolução e implementação .............................................................................7 
4. PROBLEMA DE VALOR DE CONTORNO (PVC) LINEAR .......................................11 
4.1 EQUAÇÃO DE LEGENDRE .........................................................................................11 
4.2 MÉTODO DE DIFERENÇAS FINITAS........................................................................12 
4.2.1. Resolução e implementação ...........................................................................14 
5. DISCURSÃO DOS RESULTADOS ...............................................................................17 
4. CONCLUSÃO ..................................................................................................................17 
5. BIBLIOGRAFIA...............................................................................................................18 
ANEXOS ( Algorítmos) .......................................................................................................19 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
2 
 
 
 
1. INTRODUÇÃO 
Equações diferenciais ordinárias em problemas de valor inicial (PVI) ou de contorno 
(PVC) são frequentes em modelos de diversos campos de estudo, como economia, biologia, 
mecânica dos fluidos etc. [1] Entretanto, encontrar a solução analítica para estes problemas, 
por mais que haja garantia teórica de sua existência e unicidade, pode ser demasiado 
complicado. Assim, aproximar a solução numericamente se torna conveniente. [1] 
Modelos constituem-se de abstrações, relações matemáticas que caracterizam o 
fenômeno em estudo. Visam entender ou explicar o problema em questão. Sua formulação é 
realizada seguindo passos bem definidos: selecionar as variáveis que influenciam no 
fenômeno; passar o problema para a linguagem matemática; resolução; verificar a 
consistência da solução com o problema; modificar o modelo, caso a aproximação não seja 
satisfatória; e finalmente, aplicar a solução encontrada, ou seja, ser capaz de entender o 
fenômeno por meio do modelo desenvolvido.[2] 
Este trabalho se encaixa em solucionar o modelo, que já foi fornecido em linguagem 
matemática. As soluções aproximadas primeiramente foram verificadas existentes e únicas. 
Em seguida, foram encontradas e comparadas com as soluções analíticas calculadas. 
2. FUNDAMENTAÇÃO TEÓRICA 
Dentre os diversos métodos de aproximação numérica para solução de PVC’s e 
PVI’s, foram empregados o de Runge-Kutta de 4ª ordem e Adams-Bashforth-Moulton para 
o PVI; e o de diferenças finitas para o PVC. 
Os métodos de Runge-Kutta buscam aproveitar os de Taylor, mas eliminando a 
necessidade de calcular derivadas, que é a maior desvantagem destes. Assim, visa tornar-se 
computacionalmente mais viável. Em contraste, é preciso calcular f (x , y) em diversos pontos 
[1], da seguinte forma: 
3 
 
{
 
 
 
 
 
 
 
 
𝑘1 = ℎ𝑓(𝑥𝑛0, 𝑦𝑛)
𝑘2 = ℎ𝑓 (𝑥𝑛 +
ℎ
2
, 𝑦𝑛 +
𝑘1
2
)
𝑘3 = ℎ𝑓 (𝑥𝑛 +
ℎ
2
, 𝑦𝑛 +
𝑘2
2
)
𝑘4 = ℎ𝑓(𝑥𝑛 + ℎ, 𝑦𝑛 + 𝑘3)
𝑦𝑛+1 = 𝑦𝑛 +
1
6
(𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4)
 (2.1) 
Já Adams-Bashforth-Moulton é um método de previsão-correção. Portanto, as 
primeiras aproximações de y são realizadas por um método explícito (no caso, o de Adams-
Bashforth), e posteriormente, f (x , y) é aproximada por um método implícito (Adams-
Moulton, no caso) [1] da seguinte forma: 
yn+1
(0)=yn+h/24(55 f n−59 f n−1+37 f n−2−9 f n−3) , para n=3,4,..., N 
f (n0+)1= f (xn+1, y(n0+)1) 
y(nk+)1=yn+h/24[9 f (xn+1, y(nk+−11))+19 f n−5 f n−1+ f n−2] , para k=1,2,... 
até k atingir um valor máximo pré-estabelecido ou até 
−ε<(y(nk+)1−y(nk+−11))/ y(nk+)1<ε 
 
 
 
 
4 
 
3. PROBLEMA DE VALOR INICIAL (PVI) 
𝑦′(𝑥) = 𝑦2 cos(𝑥) ∀ 𝑥 ∈ [𝜋,
3𝜋
2
] 
𝑦(𝜋) = 1 
Solução: 
Primeiramente, é necessário verificar a existência e unicidade da solução. 
Pelo Teorema de Existência e Unicidade (TEU), é preciso que as funções f (x)=y'(x) 
e ∂ f /∂ y sejam contínuas num intervalo que contenha x0 e y0 . [1] 
(i) f (x) é um produto de funções contínuas no intervalo, logo, também o é. 
(ii) f y=2 y y' cos(x) , da mesma forma produto de funções contínuas no intervalo 
de interesse, configurando-se como tal. 
Assim, estão garantidas a existência e unicidade da solução deste PVI, de forma que 
pode-se aplicar métodos, sejam analíticos ou numéricos, para encontrá-la. Resolvendo a 
equação analiticamente tem-se: 
{
 
 
 
 
 
 
 
 
 
 
 
 
 
 
𝑦′
𝑦2
= cos(𝑥)
∫
𝑦′𝑑𝑦
𝑦2
= ∫cos(𝑥) 𝑑𝑥
−
1
𝑦
= 𝑠𝑒𝑛(𝑥) + 𝑐
𝑦(𝑥) = −
1
𝑠𝑒𝑛(𝑥) + 𝑐
𝑦(𝜋) = 1 = −
1
𝑠𝑒𝑛(𝜋) + 𝑐
1 = −
1
𝑐
𝑐 = −1
 (3.1) 
Assim, a solução analítica para o PVI é 𝑦(𝑥) = −
1
𝑠𝑒𝑛(𝑥)+𝑐
 ou, de outra forma, 𝑦(𝑥) =
1
𝑠𝑒𝑛(𝑥)+𝑐
. 
5 
 
 
Figura 1. Solução analítica. 
 
 
3.1 RUNGE E KUTTA DE 4ª ORDEM 
 
3.1.1 Resolução e implementação 
 
 O algorítmo do método de Runge e Kutta baseado na referência [4] e modificado de 
acordo com o problema em questão. Inicialmente é feito apenas o cálculo do número de passos 
de acordo com o intervalo de x e o espaçamento inserido. Em seguida é obtida as soluções de 
acordo com as fórmulas para o método. Por fim é feita a plotagem da tabela e do gráfico com 
base nos valores de x e y a cada passo. 
 Foram obtidos os seguintes resultados para o espaçamento de 10-1, 10-2 e 10-3. 
6 
 
 
Figura 2. Plotagem dos valores obtidos para o passo de 0.1. 
 
Figura 3. Plotagem dos valores obtidos para o passo de 0.01. 
7 
 
 
Figura 4. Plotagem dos valores obtidos para o passo de 0.001. 
 
3.2 MÉTODO DE ADAMS-BASHFORTH-MOULTON DE 4ª ORDEM 
 
Adams-Bashforth-Moulton é um método de previsão-correção. Portanto, as primeiras 
aproximações de y são realizadas por um método explícito (no caso, o de Adams-
Bashforth), e posteriormente, f (x , y) é aproximada por um método implícito (Adams-
Moulton, no caso) [1] da seguinte forma: 
y(n
0
+
)
1=yn+h/24(55 f n−59 f n−1+37 f n−2−9 f n−3) , para n=3,4,..., N 
f(0) n+1= f (xn+1, y
(0)
(n+1)) 
y(nk+)1=yn+h/24[9 f (xn+1, y
k-1
(n+1)+19 f n−5 f n−1+ f n−2] , para k=1,2,... 
até katingir um valor máximo pré-estabelecido ou até 
−ε<(yn+1k−yk-1n+1))/ ykn+1<ε 
 
3.2.1 Resolução e implementação 
 
8 
 
Adams-Bashforth-Moulton (ABM), que é um método de previsão-correção, a ideia é 
a já mencionada: utilizar um método explícito para determinar as primeiras aproximações da 
função, pois ABM não é autoiniciável, e de posse destes valores, calcular os seguintes.Para a questão foi escolhido o método de Runge-Kutta de 4ª ordem para a 
aproximação dos primeiros valores. 
De posse destes valores, os seguintes são obtidos por previsão (método explícito de 
Adams-Bashforth) e para cada valor encontrado é feita uma correção, pelo método implícito 
de Adams-Moulton, cujas iterações devem continuar quando uma das seguintes condições 
for verdade: o erro relativo entre os valores de correção for maior que a tolerância ou a 
quantidade de iterações de correção for menor que um máximo escolhido (no caso, 10). 
Antes da aplicação do método, é necessário garantir sua convergência. As exigências 
são as seguintes: 
(i) f (x , y) e ∂ f /∂ y contínuas no intervalo; 
(ii) h escolhido de tal forma que −2<h(𝜕 f /𝜕y)<2 
(i) é verificado para a questão, como já apontado anteriormente. (ii) também é satisfeito, pois 
os h 's escolhidos são 0.1 , 0.01 , 0.001 , que satisfazem a proposição. 
9 
 
 
Figura 5. Solução Numérica com espaçamento h = 0.1, método ABM 
 
 
Figura 6. Solução Numérica com espaçamento h = 0.01, método ABM 
10 
 
 
Figura 6. Solução Numérica com espaçamento h = 0.001, método ABM 
 
11 
 
4. PROBLEMA DE VALOR DE CONTORNO (PVC) LINEAR 
 
[x(1−x)y ']'+σ(1+σ) y=0∀ x 𝜖 [0.1,0.5] 
y(0.1)=1 
y(0.5)=0 
Antes de aplicarmos qualquer método para encontrar ou aproximar a solução, da 
mesma forma que no PVC, é preciso garantir a existência e unicidade desta. O TEU 
neste caso também tem duas exigências: 
Para o problema na forma y' '= p(x) y'+q(x) y+r(x) exige-se: 
(i) p(x) limitada no intervalo. Para o 
 caso em questão tem-se p(x)=(2x−1)/(x−x2) , que 
no intervalo trabalhado é limitada por 0 
(ii) q(x) positiva em todo o intervalo. No problema trabalhado, q(x)=2/[ x(x−1)] , que é 
positiva no intervalo [−1,2] , em particular, é positiva no intervalo trabalhado. 
 
4.1 EQUAÇÃO DE LEGENDRE 
 
A equação de Legendre é dada pela expressão: 
(1 − 𝑥2)
𝑑²𝛩
𝑑𝑥²
− 2𝑥
𝑑𝛩
𝑑𝑥
+ [𝑙(𝑙 + 1) −
𝑚2
1 − 𝑥2
] 𝛩 = 0 (3.1.1) 
quando m=0. A equações aparece quando resolvida a equação de Laplace em coordeanadas 
esféricas. Esta associada à variável 𝜃, já que x = cos𝜃. Como 0 ≤ 𝜃 ≤ 𝜋, o intervalo para x é -1≤ 
x ≤ 1. 
 A equação pode ser reescrita como: 
𝑑²𝛩
𝑑𝑥²
−
2𝑥
(1 − 𝑥2)
𝑑𝛩
𝑑𝑥
+
𝑙(𝑙 + 1)
(1 − 𝑥2)
𝛩 = 0 (3.1.2) 
12 
 
 Esta equação precisa ser resolvida mediante o método de sérias. Ela tem pontos singulares 
em x = ∓1, pois estes dois valores são zeros do denominador. Isto significa que, se quisermos 
encontrar a solução em torno desses valores, precisamos utilizar algum método. No entanto, 
podemos obter a solução em torno de outro ponto, como x = 0, e, neste caso, como se trata de 
um ponto ordinário ( não é raiz do denominador, devemos utilizar o método normal de séries, 
que consiste em supor que a solução da equação diferencial é a série de potencias [5] 
𝛩(𝑥) = 𝑂(𝑥) = ∑𝑎𝑛𝑥
𝑛
∞
𝑛=0
 (3.1.3) 
 
4.2 MÉTODO DAS DIFERENÇAS FINITAS 
 
Seja a equação da forma: 
 
{
𝑦′′(𝑥) + 𝑝(𝑥)𝑦′(𝑥) + 𝑞(𝑥)𝑦(𝑥) = 𝑟(𝑥)
𝑦(𝑥0) = 𝛼
𝑦(𝑥𝑛+1) = 𝛽
 (3.2.1) 
 
onde 𝑝(𝑥) = 
𝑄(𝑥)
𝑃(𝑥)
, 𝑞(𝑥) = 
𝑅(𝑥)
𝑃(𝑥)
 e 𝑟(𝑥) = 
𝑆(𝑥)
𝑃(𝑥)
, em n partes iguais de comprimento ℎ = 
𝛽−𝛼
𝑛
, 
cada. 
Assim, 
𝑥𝑘 = 𝑥0 + 𝑘ℎ, 𝑘 = 0,1, … , (𝑛 − 1) 𝑒 𝑦𝑘 ≅ 𝑦(𝑥𝑘) = 𝑦(𝑥0 + 𝑘ℎ), 𝑘 = 0,1,2, … , 𝑛. 
Considerando a aproximação numérica para a primeira derivada utilizando a regra da 
diferença centrada: 
𝑦′(𝑥) =
𝑦𝑖+1 + 𝑦𝑖−1
ℎ
 (3.2.2) 
A segunda derivada: 
𝑦′′ =
𝑦𝑖+1 − 2𝑦𝑖 + 𝑦𝑖−1
ℎ²
 (3.2.3) 
Substituindo na equação (3.2.1), temos: 
 
(
𝑦𝑖+1 − 2𝑦𝑖 + 𝑦𝑖−1
ℎ²
) + 𝑝(𝑥𝑖)
𝑦𝑖+1 + 𝑦𝑖−1
ℎ
+ 𝑞(𝑥𝑖)𝑦𝑖 = 𝑟(𝑥𝑖) (3.2.4) 
 Para cada i, a equação discretizada fica: 
13 
 
(1 +
ℎ
2
𝑝(𝑥𝑖)) 𝑦𝑖+1 + (𝑞(𝑥𝑖)ℎ
2 − 2)𝑦𝑖 + (1 −
ℎ
2
𝑝(𝑥𝑖)) 𝑦𝑖−1 = 𝑟(𝑥𝑖)ℎ
2 (3.2.5) 
Para i=1: 
(1 +
ℎ
2
𝑝(𝑥1)) 𝑦2 + (𝑞(𝑥0)ℎ
2 − 2)𝑦1 = 𝑟(𝑥0)ℎ
2 − (1 −
ℎ
2
𝑝(𝑥1))𝛼 (3.2.6) 
Para i = n-1: 
(𝑞(𝑥𝑛−1)ℎ
2 − 2)𝑦𝑛−1 + (1 −
ℎ
2
𝑝(𝑥𝑛−1)) 𝑦𝑛−2 = 𝑟(𝑥𝑛−1)ℎ
2 − (1 +
ℎ
2
𝑝(𝑥𝑛−1))𝛽 (3.2.7) 
Assim, para determinar os pontos de 2≤ i ≤ (n-2), segue a relação descrita na equação 
(3.2.5). Logo, para resolver o sistema de equações algébricas lineares, que é um sistema de ordem 
(n-1) com matriz A tridiagonal, dada por: 
𝐴 =
[
 
 
 
 
 
 
𝑑1 𝑐1
𝑎2 𝑑2 𝑐2
𝑎3 𝑑3 𝑐3
⋱ ⋱ ⋱
𝑎𝑛−2 𝑑𝑛−2 𝑐𝑛−2
𝑎𝑛−1 𝑑𝑛−1]
 
 
 
 
 
 
 (3.2.8) 
onde 
{
 
 
 
 
𝑑𝑖 = (𝑞(𝑥𝑖)ℎ
2 − 2), 1 ≤ 𝑖 ≤ (𝑛 − 1)
𝑐𝑖 = (1 +
ℎ
2
𝑝(𝑥𝑖)) , 1 ≤ 𝑖 ≤ (𝑛 − 2)
𝑎𝑖 = (1 −
ℎ
2
𝑝(𝑥𝑖)) , 2 ≤ 𝑖 ≤ (𝑛 − 1)
 (3.2.9) 
O sistema matricial é da forma: 
A Y = B (3.2.10) 
 O vetor Y é: 
[
 
 
 
 
 
𝑦1
𝑦2
𝑦3
⋮
𝑦𝑛−2
𝑦𝑛−1]
 
 
 
 
 
 (3.2.11) 
 E por fim, o vetor B é o vetor direita da equação , na forma: 
14 
 
[
 
 
 
 
 
 
 
 
 𝑟(𝑥1)ℎ
2 − (1 −
ℎ
2
𝑝(𝑥1))𝛼
ℎ²𝑟(𝑥2)
ℎ2𝑟(𝑥3)
⋮
ℎ²𝑟(𝑥𝑛−2)
𝑟(𝑥𝑛−1)ℎ
2 − (1 +
ℎ
2
𝑝(𝑥𝑛−1))𝛽 
]
 
 
 
 
 
 
 
 
 
 (3.2.12) 
 
4.2.1. Resolução e implementação 
 
No algorítmo (Anexo III), baseado na referência [5] e modificado de acordo com o 
problema em questão, inicialmente é feita apenas o cálculo do número de passos e o 
preenchimento dos vetores x e y nos extremos. Para a avaliação dos valores durante o 
intervalo, é empregado o comando inline próprio do MATLAB R2015a. Em seguida é feito 
o preenchimento dos vetores correspondentes às três diagonais da matriz, juntamente com o 
vetor d, que é o vetor à direita da Eq.(3.2.9). Nas linhas restantes, é feita a chamada da função 
TDMAsolver (Anexo IV), algorítmo de Thomas de acordo com a referência [6], que resolve 
sistemas tridiagonais recebendo os quatros vetores do sistema. O vetor z, o qual recebe o 
resultado da função TDMAsolver. O vetor y irá auxiliar na construção da tabela e na 
plotagem do gráfico. 
O algorítmo é feito com objetivo de solucionar a Equação de Legendre deslocada: 
[𝑥(1 − 𝑥)𝑦′]′ + 𝜎(1 + 𝜎)𝑦 = 0, 𝜎 = 0.5 (3.2.1.1) 
Com base na equação (3.2.1) encontramos: 
{
 
 
 
 𝑝(𝑥) =
1 − 2𝑥
𝑥 − 𝑥²
𝑞(𝑥) =
0,75
𝑥 − 𝑥²
𝑟(𝑥) = 0
 (3.2.1.2) 
15 
 
Para o problema na forma y' '= p( x) y'+q( x) y+r (x) exige-se: 
(i) p(x) limitada no intervalo. Para o caso em questão tem-se p(x)= (2x-1)/(x-
x²), que no intervalo trabalhado é limitada por 0 
(ii) q(x) positiva em todo o intervalo. No problema trabalhado, 
q(x)=(0.75/(x(x-1)), que é positiva no intervalo [−1,1.75] , em particular, é positiva 
no intervalo trabalhado. 
 Foram obtidos os seguintes resultados para o espaçamento de 10-1, 10-2 e 10-3. 
 
Figura 8. Plotagem dos valores obtidos para o passo de 0.1. 
16 
 
 
Figura 9. Plotagem dos valores obtidos para o passo0.01. 
 
Figura 10. Plotagem dos valores obtidos para o passo 0.001. 
 
 
 
 
 
17 
 
5. DISCURSÃO DOS RESULTADOS 
Percebe-se que a solução numérica se aproxima com grande precisão da solução 
analítica, mesmo com h=0.1 , que foi o maior espaçamento empregado, tanto para Runge-
Kutta quanto para Adams-Bashforth-Moulton. Assim, todos os métodos foram convergentes 
e eficientes. Para o método ABM, os erros foram da ordem de 10−4 para h=0.1 , da ordem de 
10−8 para h=0.01 e da ordem de 10−12 para h=0.001 . Por outro lado, o tempo de execução do 
algoritmo foi perceptivelmente superior quando da diminuição do tamanho do espaçamento, 
o que pode ter sido por conta das correções. 
 
6. CONCLUSÃO 
Os métodos aplicados se mostrarm eficientes e úteis para o cálculo da aproximação 
da solução. Os erros foram pequenos, apesar de o custo computacional em termos de tempo 
de execução ter sido relativamente lento. 
 
 
 
 
 
 
 
 
 
 
 
18 
 
7. BIBLIOGRAFIA 
 
[1] RUGGIERO, M. A. G. LOPES, V. L. da R., Cálculo Numérico: aspectos teóricos e 
computacionais, 2ª edição, São Paulo: Pearson, 1997. 
 
[2] BASSANEZI, R. C. Equações Diferenciais Ordinárias: Um Curso Introdutório. 
Disponível em <http://gradmat.ufabc.edu.br/disciplinas/listas/iedo/notasdeaulas/equacoes-
diferenciais-ordinriasrodney.pdf>. Acesso em 8 de abril de 2016. 
 
[3] Código em MatLab. Disponível em 
<http://www.cfm.brown.edu/people/dobrush/am33/files/exmatlab.txt>. Acesso em 7 de 
abril de 2016. 
 
[4] Disponível em <https://marleson.wordpress.com/2010/11/03/algoritmos-em-matlab/>. 
Acesso em 04 de Abril de 2016. 
[5] MACHADO, K. D. Eletromagnetismo. Ponta Grossa, PR, Ed. TODAPALAVRA, 
2012, vol.1. 
[6] Disponível em <http://www.rc.unesp.br/igce/demac/balthazar/analise/cap9.pdf>. 
Acesso em 03 de Abril de 2016. 
[7] Disponível em < https://pt.wikipedia.org/wiki/Algoritmo_de_Thomas>. Acesso em 06 
de Abril de 2016. 
 
 
 
 
 
 
19 
 
ANEXOS (Algoritmos) 
ANEXO I 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%MÉTODO DE RUNGE-KUTTA DE ORDEM 4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%Algoritmo para calcular a solução de PVIs de ordem 1 
%Exemplo: Para calcular a solução de y'(t)= y(x)²cos(x), y(pi)=1, pi <= t <= 3pi/2. 
% f =(x,y) y(x)²cos(x) := ED a ser calculada 
% a = 0 := x inicial 
% y_inicial=1 := condição inicial 
% n = número de iterações 
% h = passo ou incremento 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% 
 
function[sol]=RK4(f,x0,y_inicial,h,n) 
 
clc 
 
f = inline ('y^2*cos(x)', 'x','y'); 
%f = ('y^2*cosx'); 
h = input('Informe o espacamento: '); 
a = input('Informe o valor inicial de x: '); 
y_inicial = input ('informe o valor inicial de y: '); 
n = ((3*pi/2)-a)/h; 
x(1)=a; 
vi= y_inicial; 
 
%Fórmulas para método de Runge-Kutta de 4a ordens 
 
for i=1:n-1, 
 
k1=h*f(a,y_inicial); 
k2=h*f(a+(1/2)*h,y_inicial+(1/2)*k1); 
k3=h*f(a+0.5*h,y_inicial+(1/2)*k2); 
k4=h*f(a+h,y_inicial+k3); 
y=y_inicial+(1/6)*(k1+2*k2+2*k3+k4); 
y_inicial=y; 
x(i+1) = x(i) + h; 
sol(i)=y_inicial; 
 
end 
 
 
% Plotagem da tabela de valores de x e y a cada passo. 
 
fprintf('\n A tabela de valores de x e y é: \n'); 
20 
 
fprintf(' k x y \n'); 
for i=1:n-1 
fprintf('%i %.3f %.4f \n', i, x(i), sol(i)); 
end 
 
sol=[vi;sol']; 
 
% Plotagem no gráfico. 
plot(x,sol,'r--'); 
legend('Solução Númerica'); 
xlabel('x'); 
ylabel('sol'); 
end 
 
 
ANEXO II 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%MÉTODO DE ADAMS-BASHFORTH-MOULTON 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function y = ABM (f, x0, y0, xn, h, er) 
 
%As entradas f, x0 e y0 compõem o PVI na forma y'(x0) = f(x,y) e y(x0) = 
y0 
%y é a aproximação da função procurada no intervalo de x0 a xn, com 
o tamanho do espaçamento h 
%er é a tolerância para a correção de Adams-Moulton 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% 
 
x(1) = x0; %primeiro valor de x, fornecido pelo usuário 
y(1) = y0; %primeiro valor de y, fornecido pelo usuário 
n = (xn-x0)/h; %quantidade de pontos nos quais y será aproximada 
 
for i = 1:3 %Runge-Kutta para determinar as primeiras 3 aproximações 
%de y, já que Adams-Bashforth-Moulton não é autoiniciável 
 
 k1 = f(x(i), y(i)); 
 k2 = f(x(i) + h/2, y(i) + h/2*k1); 
 k3 = f(x(i) + h/2, y(i) + h/2*k2); 
 k4 = f(x(i) + h, y(i) + h*k3); 
 y(i+1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4); 
 x(i+1) = x(i) + h; 
endfor 
 
for i = 4:n %método Adams-Bashforth: previsor 
 y(i+1) = y(i) + h/24*(55*f(x(i), y(i)) - 59*f(x(i-1),y(i-1)) + 
37*f(x(i-2),y(i-2)) - 9*f(x(i-3),y(i-3))); 
 x(i+1) = x(i) + h; 
 k = 0; %relativo à quantidade de iterações de correção 
 e = 1; %relativo ao erro nas iterações de correção 
 aux = y(i+1); %auxiliar para a correção 
21 
 
 
 while(k < 10 || e > er) %verificação da correção: segue enquanto as 
iterações forem 
 
 %menores que o máximo estabelecido (10) ou a tolerância for maior que 
a mÃ-nima estabelecida pelo usuário 
w = y(i) + h/24*(9*f(x(i+1),aux) + 19*f(x(i),y(i)) - 5*f(x(i-1),y(i-1)) + 
f(x(i-2),y(i-2))); %Adams-Moulton: corretor 
e = abs((w - aux)/aux); %cálculo do erro entre as iterações de 
correção 
aux = w; %auxiliar para a próxima iteração de correção 
k = k+1; %incremento das iterações 
 
 endwhile 
 
 endfor 
 
 endfunction 
 
 
ANEXO III 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%MÉTODO DE DIFERENÇAS FINITAS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [x,y] = EDO2DiferencaFinita(p, q, r, a, b, h,y_inicial, y_final) 
 
% A função calcula a EDO de Legendre deslocada de segunda ordem com valores de 
%contorno a serem inseridos utilizando o método das 
%diferenças finitas. 
% Considerando a EDO da forma: y'' + p(x)y' + q(x)y = r(x). 
% Equação de Legendre [x(1-x²)y']'+a(1+a)y = 0, p/ a = 0,5 
% As variáveis de entrada são: 
% * p = O coeficiente p(x) da equação. 
% * q = O coeficiente q(x) da equação. 
% * r = O coeficiente r(x) da equação. 
% * a = O valor inicial de x. 
% * b = O valor final de x. 
% * h = O tamanho do passo ou interações. 
% * y_inicial = O valor inicial da função y(x). 
% * y_final = O valor final da função y(x). 
% As variáveis de saída são os vetores x e y. 
% Determinação da quantidade de segmentação. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 
 
h = input('Informe o espacamento: '); 
a = input('Informe o valor inicial de x: '); 
b = input('Informe o valor final de x: '); 
22 
 
y_inicial = input ('informe o valor inicial de y: '); 
y_final = input ('informe o valor inicial de y: '); 
n = (b-a)/h; 
 
% Determinação dos vetores x e y. 
 
x(1) = a; 
x(n+1) = b; 
y(1) = y_inicial; 
y(n+1) = y_final; 
 
% Determinação das funções p, q e r pelo comando inline em funcao de x. 
 
px = inline('(1-2*x)/(x-x^2)','x'); 
qx = inline('1.75/(x-x^2)','x'); 
rx = inline('0','x'); 
 
% Determinação dos vetores diagonais da matriz tridiagonal A de coeficientes. 
% Subdiagonal 
 
a = zeros(1,n-1); 
 
% Diagonal principal 
 
b = zeros(1,n-1); 
 
% Superdiagonal 
 
c = zeros(1,n-1); 
 
% Determinação do vetor d do sistema Ay = d. 
 
d = zeros(1,n-1); 
 
% Discretização do intervalo x. 
 
for i=1:n 
x(i+1) = x(i) + h; 
end 
 
% Processo de preenchimento da subdiagonal. 
 
for j=1:n-2 
a(j+1) = 1 - (h/2)*px(x(j+2)); 
end% Processo de preenchimento da diagonal principal. 
 
23 
 
for j = 1:n-1 
b(j) = h^2*qx(x(j+1)) - 2; 
end 
 
% Processo de preenchimento da superdiagonal. 
 
for j=1:n-2 
c(j) = 1 + (h/2)*px( x(j+1)); 
end 
 
% Processo de preenchimento do vetor d. 
 
for j=1:n-1 
d(j) = rx(x(j+1))*h^2 ; 
end 
 
% Preenchimento manual dos extremos do vetor d 
 
d(1) = rx(x(2))*h^2 +((h/2)*px(x(2)) - 1)*y(1); 
d(n-1) = rx(x(n))*h^2 -((h/2)*px(x(n)) + 1)*y(n+1); 
 
% Chamada da função TDMAsolver que resolve o sistema tridiagonal. 
% Z é a solução nos pontos intermediários. 
 
z = TDMAsolver(a,b,c,d); 
 
y = [y_inicial, z(1:n-1),y_final]; 
 
% Plotagem da tabela de valores de x e y a cada passo. 
 
fprintf('\n A tabela de valores de x e y é: \n'); 
fprintf(' k x y \n'); 
for i=1:n-1 
fprintf(' %i | %.3f | %.4f \n', i, x(i), y(i)); 
end 
 
% Plotagem no gráfico. 
 
plot(x,y,'r-*'); 
legend('Solução Númerica'); 
xlabel('x'); 
ylabel('y'); 
end 
 
ANEXO IV 
 
24 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Algoritimo de Matriz Tridiagonal 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% O algorítmo de Thomas ( ou algorítmo de matriz diagonal) é um método 
% algébrico que consiste numa simplificação deeleminação gaussiana para 
% resolução de sistemas de equação tridiagonal. 
%a, b, c são os vetores colunas para a compressão da matriz tridiagonal, d é o vetor à 
direita 
% N é o número de linhas 
 
function x = TDMAsolver(a,b,c,d) 
 
N = length(d); 
 
% Modifica o primeiro coeficiente de cada linha 
c(1) = c(1) / b(1); % Risco de divisão por zero. 
d(1) = d(1) / b(1); 
 
for n = 2:1:N 
 temp = b(n) - a(n) * c(n - 1); 
 c(n) = c(n) / temp; 
 d(n) = (d(n) - a(n) * d(n - 1)) / temp; 
end 
 
% Agora voltando a substituição. 
x(N) = d(N); 
for n = (N - 1):-1:1 
 x(n) = d(n) - c(n) * x(n + 1); 
end 
end