0% acharam este documento útil (0 voto)

11 visualizações120 páginas

Cálculo Diferencial Fracionário: Fundamentos

Notas sobre cálculo fracionário, aspectos teóricos e computacionais.

Enviado por

Eduardo Castilho

Direitos autorais

Levamos muito a sério os direitos de conteúdo. Se você suspeita que este conteúdo é seu, reivindique-o aqui.

Formatos disponíveis

Baixe no formato PDF, TXT ou leia on-line no Scribd

0% acharam este documento útil (0 voto)

11 visualizações120 páginas

Cálculo Diferencial Fracionário: Fundamentos

Notas sobre cálculo fracionário, aspectos teóricos e computacionais.

Enviado por

Eduardo Castilho

Direitos autorais

Levamos muito a sério os direitos de conteúdo. Se você suspeita que este conteúdo é seu, reivindique-o aqui.

Formatos disponíveis

Baixe no formato PDF, TXT ou leia on-line no Scribd

Cálculo Diferencial e Integral Fracionário:

Aspectos Analiticos e Numéricos

Autor: J. E. Castilho
Data: Maio-2025
José Eduardo Castilho

Cálculo Diferencial e Integral Fracionário:

Aspectos Analiticos e Numéricos

Maio/2025
B
Dedico este trabalho a minha esposa Maria Cristina e filhas Alex e Laura.

C
Conteúdo

Capı́tulo 1 Breve Histórico do Cálculo Diferencial Fracionário 3

1.1 Os primeiros estágios 1695-1822 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 A integral de Abel e suas consequências. (1823-1916) . . . . . . . . . . . . . . . . 6
1.3 Cálculo Fracionário a partir de 1917 . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Capı́tulo 2 Elementos do Cálculo Clássico 12

2.1 Diferenciação e Integração . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Esquemas Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Exemplos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Capı́tulo 3 Transformadas Integrais e Funções Especiais 26

3.1 Transformadas Integrais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Função Gama de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 A Função Beta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4 Função de Mittag-Leffler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Capı́tulo 4 Cálculo Fracionário 38

4.1 Operador de Riemann-Liouville . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2 Operador de Caputo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3 Operador de Grünvald-Letnikov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Capı́tulo 5 Equação Diferencial Ordinária Fracionária 54

5.1 Equações do tipo Riemann-Liouville e do tipo Caputo . . . . . . . . . . . . . . . . . 54
5.2 Métodos de Diferenças Fracionárias Regressivas . . . . . . . . . . . . . . . . . . . 64
5.3 Método do tipo Adams-Bashforth-Moulton . . . . . . . . . . . . . . . . . . . . . . . 69
5.4 Exemplos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Capı́tulo 6 Equação Diferencial Parcial Fracionáriado tipo Riesz 76

6.1 Equação de Advecção-Difusão Fracionária do tipo de Riesz . . . . . . . . . . . . . . 76
6.2 Esquemas Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.3 Exemplos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Capı́tulo A Esquema Adams-Bashforth 88

Capı́tulo B Esquema Adams-Moulton 90

Capı́tulo C Esquema de Adams-Bashforth-Moulton 92

CONTEÚDO

Capı́tulo D Esquema de aproximação L1-L2 95

Capı́tulo E Esquema de aproximação Grünwald deslocado 102

Capı́tulo F Esquema de aproximação por Transformação Matricial 106

Bibliografia 110

ii
Introdução
Muitos autores consideram que o surgimento do cálculo diferencial fracionário decorreu da troca
de cartas entre L’Hopital e Leibniz em 1695. Na ocasião L’Hopital questionava Leibniz sobre o
significado de derivadas não inteiras, especialmente o caso d1/2 xy ( Notação da época para a derivada
√
de y(x) = x ). Leibniz assegurava que a resposta para este caso seria d1/2 x = x dx : x. Dizia que
“um dia consequências muito úteis serão tiradas desse paradoxo, já que existem pequenos paradoxos
sem utilidade” [Oldhann74]. O cálculo de ordem não inteira tem a intensão de generalizar o cálculo
classico de ordem inteira, proposto por Newton e Leibniz, no século XVII. Popularmente é chamado
de Cáculo Fracional ou Cálculo Fracionário, mas nos dias hoje tem-se uma visão mais atualizada,
que considera a integração e diferenciação de ordem arbitrária, diferentemente da visão restrita que o
enxerga apenas como cálculo de frações. Neste texto será adotado o termo popular Cálculo Fracionário,
já que a abordagem numérica será feita em ordem fracionária.
No desenvolvimento do cálculo clássico, surge classes de funções especiais nas soluções de
equações diferenciais ordinárias e parciais. No Cálculo Fracionário, também temos diversas funções
especiais, tais como a Função Gama, que generaliza o conceito de fatorial, a Função de Mittag-Leffler
de um parâmetro, como a generalização da função exponencial.
Existe mais de uma formulação possı́vel para o Cálculo Fracionário, sendo cada uma mais
adequada a um certo contexto fı́sico que outro, mas todas podem ser pensadas como operadores que
representam funções da memória sobre o processo de alguns sinais de sistemas fı́sicos. A ideia de
que fenômenos fı́sicos, como difusividade anômala ou processos ondulatórios, pode ser descrito com
modelos diferenciais fracionários levantando, pelo menos, as seguintes três questões fundamentais: Os
modelos matemáticos com derivadas fracionárias de espaço e/ou tempo são consistentes com as leis
fundamentais e simetrias bem conhecidas da natureza? Como a ordem fracionária de diferenciação
pode ser observada experimentalmente ou como uma derivada fracionária emerge de micro-modelos
escópicos? Uma vez que um modelo de cálculo fracionário esteja disponı́vel, como um modelo
fracionário pode ordenar que a equação seja resolvida (exatamente ou aproximadamente)?
O Cálculo Fracionário se manteve relativamente desconhecido por se limitar a poucas contribuições
e aplicações nas ciências aplicadas, como nas engenharias. O tema começou a despertar interesse da
comunidade cientı́fica, em 1974, com a realização da primeira conferência internacional sobre Cálculo
Fracionário, promovida por Bertram Ross em New Havem, USA. Na década de 90 ocorre uma ex-
plosão de publicações de livros e artigos sobre o tema, consolidando o Cálculo Fracionário como
uma importante ferramenta na análise de fenômenos, demonstrando ter propriedades e caracterı́sticas
superiores ao cálculo integro-deferencial clássico.
Estas notas procura abordar vários aspectos do Cálculo Fracionário estruturadas das seguintes
formas:
O Capı́tulo 1 aborda uma breve revisão histórica da teoria do Cálculo Fracionário e suas
aplicações; No Capı́tulo 2 são apresentados alguns resultados do cálculo clássico de ordem inteira, no
contexto analı́tico e numérico; Capı́tulo 3 trata de alguns resultados de operadores integrais e funções
CONTEÚDO

especiais; O Capı́tulo 4 aborda diferentes formulações do Cálculo Fracionário e suas aplicações. No

Capı́tulo 5 é apresentado alguns exemplos de equações diferenciais ordinárias fracionárias no contexto
analı́tico e numérico. A Equação de Advecção-Difusão Fracionária é analisada no Capı́tulo 6.

2
Capı́tulo 1 Breve Histórico do Cálculo Diferencial
Fracionário
O Cálculo Fracionário surgiu quase simultaneamente ao cálculo clássico, mas ganhou reconhe-
cimento como área de interesse da comunidade cientı́fica, principalmente nos anos 90. Trabalhos
recentes tem considerado a integração e diferenciação de ordem arbitrária, não se restringindo apenas
as ordens fracionárias.
Nos últimos anos, diversos artigos ([44], [43]) e livros ([30], [45]) têm abordado vários aspectos
da história e desenvolvimento do cálculo fracionário. Nas citações históricas iremos manter a notação
que era usada na época.
Com o objetivo de fornecer uma visão geral do seu desenvolvimento histórico, apresentamos uma
retrospecção dividida em três seções, abordando os diferentes estágios do desenvolvimento do cálculo
fracionário.

1.1 Os primeiros estágios 1695-1822

A origem do cálculo fracionário pode ser atribuı́da a Leibniz. Da sua notação para a n-ésima
dn y
derivada de uma função, dxn com n sendo um número inteiro, surge a ideia de considerar os casos em

que n não seja um número inteiro. Embora não haja artigos matemáticos de Leibniz sobre o assunto,
há uma série de cartas enviadas a L’Hopital, nas quais procura dar um significado a derivadas não
inteiras, especialmente o caso d1/2 xy. Em sua carta Leibniz escreve: [43]

“Você pode ver aqui, senhor, que se pode expressar um termo como d1/2 xy ou d1:2 xy por
uma série infinita, ainda que pareça distante da geometria, que costuma considerar apenas
as diferenças de expoentes inteiros positivos ou os negativos com respeito a somas, mas
não aquelas cujos expoentes são fracionários. É verdade que é ainda mostrar o que é esta
série para d1:2 x; mas não só isso pode ser explicado de certa forma. Como as ordenadas
x são expressas em uma série geométrica, tal que escolhendo uma constante dβ segue
que dx = xdβ : a, ou (se alguém escolheu a como unidade) dx = xdβ, significando
2 3 e
que ddx seria xdβ , e d3 x seria igual a xdβ etc. e de x seria igual a xdβ . E assim
o expoente diferencial foi alterado e pela substituição de dβ por dx : x, resultando em
e √
de x = dx : x x. Portanto segue-se que d1/2 x será igual a x dx : x. Parece que um
dia, consequências muitos uteis serão tiradas desse paradoxo, uma vez que há pequenos
paradoxos sem utilidade.”

Leibniz reconhece a questão levantada por L’Hospital, considera suas implicações em termos
matemáticos. O tema de derivada não inteira é retornada na correspondência entre Johann Bernoulli
e Leibniz. Em dezembro de 1695, Bernoulli reitera o tema das derivadas “fracionária ou irracional”.
Leibniz aborda esse problema em uma carta de volta a Johann Bernoulli escrita no mesmo mês, onde
1.1 Os primeiros estágios 1695-1822

menciona derivadas de “ordem geral” e repete essa questão em mais detalhes em comparação com a
carta ao L’Hospital:

“O que você procura nas derivadas, cujo expoente é fracionário ou irracional, já con-
siderei em carta ao L’Hospital. Eu acrescentei isso, as derivadas estruturadas podem
ser comparadas com equações diferenciais. Seja por exemplo d1:2 x a derivada proposta.
Permita que x seja uma progressão geométrica. Deixe a constante derivada assumida ser
dh, e temos xdh : a = dx, e obtemos d2 x = dxdh : a = xdhdh : aa e similarmente
d3 x = xdh3 : a3 , e mais geral de x = xdhe : ae , bem como d1:2 x = xdh1:2 : a1:2 , ou
√
d1:2 x = x dh : a (aqui 1 : 2 é igual a 1/2; e dh : a igual a dh
a
). Assim, você pode ver que
tais derivadas estruturadas só têm significado extraindo ou exponenciando as derivadas
ordinárias. Eu acho que isso é memorável e você não será ingrato por isso. Você viu que
as derivadas ordinárias podem ser expressas compondo séries geométricas infinitas e eu
suspeito que isso será verdade para os casos que não são reais.”

O tópico das derivadas de ordem não inteira não termina com a morte de Leibniz em 1716. Em
1783, Leonhard Euler trabalhou em progressões de números e pela primeira vez introduziu a função
Gama como generalização de fatoriais [9]:

“27. Para encerrar esta discussão, deixe-me acrescentar algo que certamente é mais
curioso do que útil. Sabe-se que dn x denota a diferencial de x de ordem n, e se p denota
qualquer função de x e dx é considerado constante, então dn p é homogêneo com dxn ; mas
sempre que n é um número inteiro positivo, a razão de dn p com dxn podem ser expressos
algebricamente; por exemplo, se n = 2 e p = x3 , então d2 (x3 ) está para dx2 como 6x
para 1. Agora perguntamos, se n é um número fracionário, qual o valor dessa proporção
deve ser. É fácil entender a dificuldade neste caso; pois se n é um número inteiro positivo,
dn pode ser encontrado por diferenciação contı́nua; tal abordagem não está disponı́vel
se n for um número fracionário. Mas ainda assim será possı́vel desembaraçar a questão
usando a interpolação em progressões, que discuti neste ensaio.”

Indiretamente Lagrange, contribuiu para o cálculo fracionário, quando desenvolveu a lei dos
expoentes para operadores diferenciais de ordem inteira [16]
Pouco mais de meio século após a morte de Leibniz, a obra de J.L. Lagrange indiretamente
contribuiu para o campo do cálculo fracionário. Em 1772 Lagrange desenvolveu a lei dos expoentes
para operadores diferenciais de ordem inteira
dm dn dm+n
= m+n , m, n ∈ N
dxm dxn dx
Este resultado pode ser generalizado para escolhas arbitrárias de n, m ∈ C sob certas condições, que
se tornou evidente muito mais tarde na história.
A primeira definição detalhada de uma derivada fracionária é atribuı́da a P. S. Laplace [17], em
que se define uma derivada fracionária para funções representáveis por uma integral:

4
1.1 Os primeiros estágios 1695-1822

Para esclarecer isso, vamos considerar a equação

Z
yx = T dtt−x
onde T é uma função de t e a integral é calculada num dado intervalo. Seja x uma
“variente” de α Z
−x 1
∆yx = T dtt −1
tα
ou generalizando
Z i
i −x 1
∆ yx = T dtt −1
tα
; sendo i negativo, o ∆ na integral altera o sinal da integral. Se α é infinitamente pequeno
e igual a dx obtém-se 1/t = 1 + dx log(1/t). Assim, observando que, ∆i yx obtém-se
i
di y x
Z
−x 1
= T dtt log
dxi t
Da mesma forma que se obtém adotando a notação de no 2

Z i
i −x b q
∇ yx = T dtt a + + ··· + n
t t

Assim, a mesma análise, que resulta em derivação bem-sucedida das variáveis das funções
R
geradoras, as funções sob , a integral definida, expressando essas derivadas. O sı́mbolo
∇i expressa estritamente apenas um conjunto i, contendo operações consecutivas. A
análise da geração das funções reduzem essas operações ao aumento de um polinômio em
seus expoentes; e o exame das integrais definidas produz diretamente a expressão ∇i yx ,
mesmo no caso de i ser um número fracionário.

No mesmo livro, Laplace generaliza seus resultados para outro conjunto de funções:

Se a função ys , dependente em s, pode ser representado por uma integral da forma

R s
x φdx, as diferenças infinitesimais de ordem n são dadas por no 21 como
dn y s
Z
= xs φdx(log x)n ,
dsn
Z
n
∆ ys = xs φdx(x − 1)n .
R
Se ao invés de representar a função de s como integral xs φdx, for usado a representação
R −sx
c φdx, obtém-se
dn ys
Z
n
= (−1) n
xn φdxc−sx ,
ds
Z
n
∆ ys = φdxc−sx (c−x − 1)n .

Para obter integrais negativas, as quais são assumidas ser finita e infinitesimal, é suficiente
escolher n negativo nas fórmulas acima. Pode-se verificar que a fórmula também é válida
para n geral, mesmo que n seja fracionário.

5
1.2 A integral de Abel e suas consequências. (1823-1916)

S. F. Lacroix, no seu livro “Traité du calcul différentiel et du calcult intégral”, 1819, elabora uma
generalização de uma derivada de ordem inteira para ordem fracionária. Sendo n-ésima derivada da
função y(x) = xm , m ∈ N dada por
dn n!
n
y(x) = xm−n , m ≥ n.
dx (m − n)!
Lacroix faz a generalização do fatorial, usando a função Gamma de Euler, para formular a derivada
não inteira, sendo n ∈ R, segue:
dn Γ (m + 1)
n
y(x) = xm−n .
dx Γ (m − n + 1)
Para o caso particular y(x) = x e n = 1/2, Lacroix apresenta a representação
√
d1/2 Γ (2) 1/2 2 x
y(x) = x = √ .
dx1/2 Γ (3/2) π
A relevância matemática do trabalho de Lacroix reside no fato de que seu resultado é a mesma
representação da atual derivada fracionário a do tipo Riemann-Liouville.
Uma definição mais geral foi apresentada por Fourier, em 1822, no seu trabalho “Théorie analy-
tique de la chaleur”. Ele observou que a expressão:

π iπ iπ
cos r + i = cos r cos − senrsen ,
2 2 2
que resulta em −senr, − cos r, senr, cos r, . . ., quando i = 1, 2, 3, 4, . . .. Desta forma função
Z +∞ Z +∞
1
f (x) = f (α)dα cos(p(x − α))dp
2π −∞ −∞
tem sua i-ésima derivada dada por:
Z +∞ Z +∞
di 1 i
π
f (x) = f (α)dα p cos p(x − α) + i dp
dxi 2π −∞ −∞ 2
Ao considerar que i pode assumir qualquer valor positivo ou negativo, Fourier fornece uma definição
para a derivada não inteira. Os trabalhos de Fourier marca o final do primeiro estágio, em que se
deram formulações para a derivada fracionária, mas sem apresentar aplicações da mesma.

1.2 A integral de Abel e suas consequências. (1823-1916)

O Cálculo Fracionário surgiu junto com o cálculo clássico, mas não foi até 1823, que as operações
fracionárias fossem usadas para resolver um problema fı́sico especı́fico. Foi Niels Henrik Abel quem
usou esta nova ferramenta matemática para resolver uma equação integral surgido no problema da
tautócrona [30].
O problema consiste em determinar uma curva no plano (x, y), onde uma partı́cula, sob a ação da
da gravidade, deve percorrer a curva no mesmo tempo. Independente da sua posição inicial (x0 , y0 )
na curva. Usando o conceito da fı́sica, que a energia potencial perdida durante a descida da partı́cula
é igual à energia cinética que a partı́cula ganha, obtemos a seguinte equação diferencial
2
1 dΛ
m = mg(y0 − y),
2 dt

6
1.2 A integral de Abel e suas consequências. (1823-1916)

Figura 1.1: Curva tautócrona com T = 10s.

onde m é a massa da partı́cula, Λ é a distância que a partı́cula percorre sobre a curva, desde o ponto
inicial e g a força gravitacional. Separando as variáveis tempo e espaço temos
dΛ p
−√ = 2gdt
y0 − y
e integrando sobre a variação de tempo t = 0 a t = T , segue
p Z y0
2gT = (y0 − y)−1/2 dΛ.
0
Dado que o tempo que a partı́cula necessita para atingir o ponto mais baixo da curva é constante,
podemos considerar que o lado esquerdo da equação acima será uma constante, que chamaremos de
k. Definindo a distância percorrida como uma função da altura (posição em y) por Λ = F (y) seque
que dλ/dt = F ′ (y). Fazendo a mudança de variáveis y0 por x e y por t e denotando F ′ = f a integral
acima se torna Z x
k= (x − t)−1/2 f (t)dt,
0

onde f é a função a ser determinada. Ao multiplicar ambos os lados por 1/Γ (1/2) Abel obtém
Z x
k 1 −1/2 d−1/2
= (x − t) f (t)dt = −1/2 f (x).
Γ (1/2) Γ (1/2) 0 dx
Usando a derivada inversa a esquerda, ou seja
d1/2 d−1/2 d0
f (x) = f (x) = f (x)
dx1/2 dx−1/2 dx0
obtemos a solução do problema da tautócrona
1 d1/2 k
f (x) = 1/2
k= √ ,
Γ (1/2) dx π x
onde a última igualdade se deve a derivada de ordem 1/2 da constante k. Na Figura (1.1) temos a
curva, tomando o valor de T = 10s. Samko, em seu livro [45] destaca que é importante observar que
Abel, não só resolve o problema de tautócrona, mas fornece uma solução para uma equação integral
mais geral, Z x
f (t)
k(x) = α
dt, x > a, 0 < α < 1.
a (x − t)

Em uma série de oito artigos, J. Liouville aborda o Cálculo Fracionário. Em [22] desenvolveu
duas definições diferentes de derivadas fracionárias. A primeira é aplicada em funções f (x) que

7
1.2 A integral de Abel e suas consequências. (1823-1916)

podem ser expressas em séries da forma

∞
X
f (x) = ck eak x ,
k=0
em que, partindo da derivadas de ordem inteiras Dn eax = an eax , substitui-se n ∈ N por α ∈ C,
obtendo-se: ∞
X
α
D f (x) = ck aαk eak x .
k=0

A escolha de α se restringe aos casos em que a série converge. A segunda definição, feita para funções
do tipo f (x) = 1/xa , é dada por:
(−1)α Γ (a + α) −a−α
Dα x−a = x ,
Γ (a)
onde a, α ∈ C. Liouville aplica essas definições em problemas geométricos, fı́sicos e mecânicos [22].
Em 1847, G. F. B. Riemann procurou generalizar a série de Taylor, usando a integral fracionária
de ordem α da função f (x) definida por:
Z x
α 1
D f (x) = (x − t)(α−1) f (t)dt + ψ(x).
Γ (α) c
A função ψ(x) foi introduzida pela ambiguidade do extremo inferior c, que corresponde a função
complementar definida no trabalho de Liouville [21]. Quando c = 0 e sem a função complementar, a
definição acima é conhecida como Integral Fracionária de Riemann-Liouville, uma das definições mais
usadas nos dias de hoje e será abordada no Capı́tulo 4. A representação, sem a função complementar
não é uma trivialidade, como observado por A. Cayley [4]. A resolução do problema da função
complementar se deve a contribuições de várias pessoas. Em 1869 N. Y. Sonin [46] usa a fórmula
da integral de Cauchy como ponto de partida para obter a diferenciação com ı́ndice arbitrário. A. V.
Letnikov estendeu a ideia de Sonin pouco tempo depois, em 1872 em seu artigo [19]. Ambos tentaram
definir derivadas fracionárias utilizando um contorno fechado. Começando com a fórmula integral de
Cauchy para derivadas de ordem inteira, dada por
Z
(n) n! f (t)
f (z) = dt,
2πi C (t − z)n+1
e a generalização para o caso de ı́ndice arbitrário é obitido trocando o fatorial pela função Gama de
Euler, onde α! = Γ (α + 1). A extensão direta para valores não inteiros resulta no problema de que o
integrando contém um ponto de ramificação, como demostrado por Ross [44]. Foi no ano de 1884 que
H. Laurent usou um contorno dado como um circuito aberto, em vez de um circuito fechado usado
por Sonic e Letnikov [18], obtendo a definição integral fracionária de Riemann-Liouville.
Z x
−α 1
c Dx f (x) = (x − t)α−1 f (t)dt, Re(α) > 0.
Γ (α) c
Quase simultaneamente ao trabalho de Sonin, Grünwald [11] e Letnikov [20] forneceram a base
para outra definição de derivada fracionária que também é frequentemente usada nos dias de hoje.
Usando diferenças de ordem fracionária, Grünwald e Letnikov obtiveram como definição da derivada
fracionária
GL α (∆αh f )(x)
Dx f (x) = lim (1.1)
h→0 hα

8
1.3 Cálculo Fracionário a partir de 1917

onde, as diferenças finitas (∆αh f )(x) são definidas como

∞
k α
X
α
(∆h f )(x) = (−1) f (x − kh), α > 0
k=0
k
α

e k é a generalização dos coeficientes binomiais, onde os fatoriais são substituı́dos pela função
Gama de Euler. Nos dias de hoje, a definição de Grünwald-Letnikov da derivada fracionária é usada
nos esquemas numéricos, nos quais se usam a forma truncada da série (1.1).
No final do século XIX, O. Heaviside [13] apresenta o cálculo operacional. Baseado na ideia
de que a operação diferencial d/dx é substituı́da por uma letra p e tratado como uma constante na
solução de equações diferenciais. Heaviside usou potências arbitrárias de p, principalmente p1/2 , para
obter soluções em problemas de engenharia. Segundo Ross [44], Heaviside interpretou p1/2 aplicado
em 1, isto é D1/2 1, como sendo 1/(πt)1/2 . Desde que f (x) = 1 é uma função da classe de Riemann,
o operador de Heaviside pode ser interpretado como o operador de Riemann 0 Dxα e desta forma,
mostrando uma aplicação do Cálculo Fracionário.
Neste ponto encerramos a segunda fase, que se caracteriza pela aplicações do cálculo fracionário
a problemas fora da matemática pura, iniciado com o problema da tautocrona de Abel.

1.3 Cálculo Fracionário a partir de 1917

Iniciamos este perı́odo fazendo referência ao trabalho de H. Weyl [51]. Weyl considera a
transformada de Fourier da função periódica ϕ(x), dada por
+∞ Z 2π
X
ikx 1
ϕk e , ϕk = eikx ϕ(x)dx
−∞
2π 0

a fim de definir a integração fracionária adequada para essas funções, dada por
Z 2π
α 1 α
I± ϕ(x) = ψ± (x − t)ϕ(t)dt, (1.2)
2π 0
α
com algumas funções especiais ψ± (x). Além disso Weyl demostrou que essas integrais fracionárias
podem ser escritas como
Z x
α 1
I+ ϕ(x) = (x − 1)α−1 ϕ(t)dt, (1.3)
Γ (α) −∞
Z +∞
α 1
I− ϕ(x) = (x − 1)α−1 ϕ(t)dt, (1.4)
Γ (α) x
dado que as integrais sejam convergentes. Weyl garante a convergência para funções ϕ(x) periódicas,
cujo o coeficiente de Fourier ϕ0 = 0. Conforme observado por Samko et al. [10] as definições atuais
de integrais de Weyl geralmente não mencionam esse detalhe. Como a definição de Riemann-Liouville
já considerava o intervalo de integração infinito, nomear as equações (1.3) e (1.4) como integral de
Weyl é um erro histórico.
Em 1918, O’Shaughnessy [33] apresenta uma solução explicita para a equação diferencial fra-
cionária Dα y = y/x. O problema não foi rigorosamente definido, desde que não se fez menção que
tipo de operador diferencial fracionário estava sendo aplicado. Tanto que, em 1919, Post [36] propõe
uma solução completamente diferente e correta.

9
1.3 Cálculo Fracionário a partir de 1917

Em 1927, Marchaud [28] desenvolve uma versão da integral de Grünwald-Letnikov, definida em

(1.1), usando Z +∞
M α (∆ℓt f )(x)
D f (x) = c dt, α > 0, (1.5)
0 t1+α
como derivada francionária de uma função f . O termo (∆ℓt f (x) é uma diferença finita de ordem
ℓ > α e c é uma constante de normalização. Sobre certas condições, a definição (1.5) coincide com
as versões da derivada fracionário de Riemann-Liouville e Grünwald-Letnikov.
No ano de 1931, Watanabe [50] desenvolve a versão da fórmula de Leibniz para a derivada
fracionário de Riemann-Liouville, dada Portanto
+∞
α
X α
D (f g) = Dα−β−k f Dβ+k g, β ∈ R (1.6)
k=−∞
k+β
onde f e g são funções analı́ticas.
Os artigos de M. Riesz [38–41] são centrados na integral
Z +∞
R α 1 ϕ(t)
I ϕ= dt, Re(α) > 0, α ∈ /N (1.7)
2Γ (α) cos(απ/2) −∞ |t − x|1−α
conhecida como potencial de Riesz. Esta integral está relacionada com as integrais de Weyl (1.3) e
(1.4 ) e consequentemente com a integral fracionária de Riemann-Liouville, onde
I = (I+α + I−α )(2 cos(απ/2))−1
R α

Uma modificação da integral fracionária de Riemann-Liouville foi apresentada nos trabalhos de

Erdélyi e Kober [7, 8], onde
2x2(α+η) x 2
Z
α,η
Ix ϕ(x) = (x − t2 )α−1 t2η+1 ϕ(t)dt, (1.8)
Γ (α) 0
2x2η ∞ 2
Z
α,η
I∞ ϕ(x) = (t − x2 )α−1 t1−2α−2η ϕ(t)dt, (1.9)
Γ (α) 0
(1.10)
que se tornaram úteis em várias aplicações. Enquanto essas ideias estão fortemente ligadas à
√
diferenciação fracionária das funções x2 e x, já feito por Liouville 1832 em [22], o fato de Erdélyi
e Kober usarem a transformada de Mellin por seus resultados é digna de nota.
Durante o Século XX, um grande número de resultados adicionais foram apresentados. Dentre
eles um merece destaque por ser amplamente utilizado. M Caputo [2]. Dada uma função f que seja
n − 1 absolutamente derivável e contı́nua, Caputo define a derivada fracionária por
Z t n
C α 1 n−α−1 d
D∗ f (x) = (t − s) f (s)ds. (1.11)
Γ (n − α) 0 ds
Essa formulação é conhecida como a derivada fracionário de Caputo. Está fortemente conectada com
a derivada de Riemann-Liouville e hoje é frequentemente utilizada nas aplicações. Isto se deve pelo
fato de que, usando a derivada de Caputo, pode-se especificar as condições iniciais de uma equação
diferencial fracionário, na forma clássica, isto é y (k) (0) = bk com k = 0, 1, 2, . . . , n − 1, fato que não
ocorre com a derivada de Riemann-Liouville.
Na segunda metade do século XX, o campo do Cálculo Fracionário havia crescido a tal ponto
que, em 1974, a primeira conferência tratou exclusivamente da teoria e aplicações na área foi realizado

10
1.3 Cálculo Fracionário a partir de 1917

em New Haven [42].

No mesmo ano Oldham e Spanier publicam o primeiro livro sobre cálculo fracionário [Oldhann74].
Desde então diversos livros tem sido publicado, sendo um dos mais populares o de Miller e Ross [30].
Desde seu surgimento, com L’Hospital questionando Leibniz, até a ampla aplicação em diversos
campos cientı́ficos, o Cálculo Fracionário tem percorrido um longo caminho. Mesmo sendo tão antigo
quanto ao próprio cálculo clássico, somente nas últimas decadas encontrou aplicabilidade em modelos
que descrevem problemas complexos da vida real.

11
Capı́tulo 2 Elementos do Cálculo Clássico
O Cálculo é uma campo da matemática bem conhecido e tem seus resultados aplicados em
diversas áreas do conhecimento. Existe diversos livros que fazem abordagem analı́tica ou numérica
sobre o cálculo diferencial e integral, análise, equações diferenciais ordinárias, equações diferenciais
parciais, equações integrais, etc. Todos esses temas podem ser abordados no contexto do Cálculo
Fracionário . Desta forma, resultados relevantes do cálculo clássico são abordados neste capı́tulo.
Vamos considerar os aspectos analı́ticos e numéricos para uso nos próximos capı́tulos.

2.1 Diferenciação e Integração

O objetivo do Cálculo Fracionário é generalizar os conceitos do cálculo clássico e manter a
relação entre os conceitos. Apresentamos alguns resultados do cálculo clássico que serão necessários
para os próximos capı́tulos. Iniciamos essa seção apresentando o teorema que faz a conexão entre o
cálculo diferencial e o cálculo integral, o Teorema Fundamental do Cálculo [48].
Teorema 2.1 (Teorema Fundamental do Cálculo)
Seja a função f : [a, b] → R contı́nua em [a, b] e seja F : [a, b] → R definida por
Z x
F (x) = f (t)dt
a
Então F é diferenciável e
d
F (x) = f (x)
dx ♡

Usaremos as seguintes notações, tanto no caso clássico e fracionário.

Observação
1. Denotamos o operador que mapeia a função diferenciável a sua derivada por:
d
Df (x) := f ′ (x) = f (x).
dx
2. Seja a função f Rieamann integrável no intervalo compacto [a, b]. Denotamos o operador que
mapeia f a sua função primitiva por
Z x
Ia f (x) := f (t)dt
a
para a ≤ x ≤ b.
3. Para n ∈ N, segue que Dn := DDn−1 e Ian := Ia Ian−1
Um resultado importante na formulação do Cálculo Fracionário é a formulação da aplicação
recursiva de Ia para uma fórmula explicita para o operador Ian , n ∈ N, que pode ser vista como uma
convolução de Laplace entre f(x) e a função de Gel’fand-Shilov de ordem n, definida como
2.1 Diferenciação e Integração

Definição 2.1 (Função de Gel’fand-Shilov)

Sejam n ∈ N. Defini-se a função de Gel’fand-Shilov como
 n−1
t
 (n−1)! se t ≥ 0

ϕn (t) :=

0 se t < 0.

♣

Lema 2.1
Seja f uma função Rieamann integravel em [a, b]. Para x ∈ [a, b] e n ∈ N temos
Z x
n 1
Ia f (x) = ϕn (x) ∗ f (x) = (x − t)n−1 f (t)dt.
(n − 1)! a
onde ∗ representa a convolução de Laplace.
♡

Demonstração Demonstra-se o Lema por indução no parâmetro n. Para n = 1 temos que, pela
definição (2.1) e a convolução de Laplace, que
Z x
(x − t)1−1
Ia f (x) = ϕ1 (x) ∗ f (x) = f (t)dt
a (1 − 1)!
Supondo que seja válido para n, seque que
Ian+1 f (x) = Ia [Ian f ](x)
= Ia (ϕn ∗ f )(x)
Z x
= ϕn (t) ∗ f (t)dt
a
Z xZ u
(u − t)n−1
= f (t)dtdu.
a a (n − 1)!
Z xZ x
(u − t)n−1
= f (t)dudt.
a u (n − 1)!
Z x
1
= (x − t)n f (t)dt
(n)! a
onde se aplica o Teorema de Goursart [23], para a troca da ordem de integração
Através desse Lema, podemos fazer uma relação entre o operador diferencial de ordem m com o
operador integral. Segue que DIa f = f , o que implica que Dn Ian f = f . Com isto temos o seguinte
Lema.
Lema 2.2
Seja m, n ∈ N com m > n, e seja f uma função continuamente derivável em [a, b]. Então,
Dn f = Dm Iam−n f
♡

Demonstração Sendo m > n seque que m − n = k > 0 e m = n + k, assim segue

Dm Iam−n f = Dn+k Iak f = Dn Dk Iak f = Dn f

Antes de prosseguir com as propriedades da integração e diferenciação, vamos introduzir alguns

espaços de funções nos quais os próximos resultados serão formulados.

13
2.1 Diferenciação e Integração

Definição 2.2
Seja 0 < µ ≤ 1, k ∈ N e p ≥ 1. Definimos
Z b
p
Lp [a, b] := f : [a, b] → R; f é mensurável em [a, b] e |f (x)| dx < ∞
a
L∞ [a, b] := {f : [a, b] → R; f é mensurável e essencialmente limitada em [a, b]}
Hµ [a, b] := {f : [a, b] → R; ∃c > 0, ∀x, y ∈ [a, b] : |f (x) − f (y)| ≤ c|x − y|µ }
C k [a, b] := {f : [a, b] → R; f tem a k-ézima derivada contı́nua }
C[a, b] := C 0 [a, b]
H0 [a, b] := C[a, b]
♣

Os espaços Lp [a, b], com 1 ≤ p ≤ ∞ são chamados de Espaços de Lebesque, enquanto Hµ [a, b]
são chamados de Espaços de Hölder ou Lipschitz de ordem µ. Um resultado que será necessário
para os próximos é a versão do Teorema Fundamental do Cálculo no Espaço de Lebesque, dado pelo
Teorema:
Teorema 2.2 (Teorema Fundamental no Espaço de Lebesgue )
Seja f ∈ L1 [a, b]. Então Ia f é diferenciável em quase todo os pontos em [a, b], e DIa f = f
também vale em quase todo o intervalo [a, b].
♡

Uma versão da demostração deste Teorema pode ser encontrada no trabalho de R. L. Pouso [37].
No cálculo clássico, são bem conhecidas as seguintes regras de derivação. Sendo f, g, ∈ C[a, b]
temos:
Regra do Produto:
D(f · g)(x) = (Df )(x) · g(x) + f (x) · (Dg)(x).

Regra do Quociente: Sendo g(x) ̸= 0

f (Df )(x)g(x) − f (x)(Dg)(x)
D (x) = .
g [g(x)]2
Regra da Cadeia:
d
Dg(f (x)) = g(f (x))Df (x).
df
O uso repetido da regra do produto fornece a Fórmula de Leibniz,
n
n
X n
D (f · g)(x) = Dk f (x) · Dn−k g(x)
k=0
k
com n ∈ N e f, g ∈ C n [a, b].
Existe uma generalização para a regra da cadeia, conhecida como Fórmula de Faà di Bruno [47].
Se f e g tem derivadas até ordem n, então
n (1) k1 (2) k2 (n) kn
n
X X n! (k) f (x) f (x) f (x)
D g(f (x)) = g (f (x))· ··· ,
k 1 !k2 ! · · · kn ! 1! 2! n!
k=1 k +k +···+k =k
1 2 n
k1 +2k2 +···nkn =n
onde as condições sobre k1 , k2 , . . . kn determinam o número de partições de {1, 2, . . . n}.

14
2.1 Diferenciação e Integração

Iremos abordar uma versão não clássica do Teorema de [Link] isto necessitamos da definição
de outro espaço de funções:
Definição 2.3
Denotamos por An ou An [a, b] o conjunto de funções que apresentam (n − 1) derivadas
absolutamente contı́nuas, isto é, funções f para qual existe, em quase toda parte, uma função
g ∈ L1 [a, b] tal que Z x
f (n−1) (x) = f (n−1) (a) + g(t)dt.
a

Neste caso, a função g é chamada de n-ésima derivada (de forma generalizada) de f , simplifi-
cando g = f (n) . ♣

Teorema 2.3
Para m ∈ N as seguintes afirmações são equivalentes:
1. f ∈ Am [a, b]
2. Para todo x, y ∈ [a, b]
m−1
X (x − y)k
f (x) = Dk f (y) + Iym Dm f (x).
k=0
k!
♡

A forma clássica do Teorema de Taylor é obtida com y = a e usando o fato de que 1. ⇒ 2.. Se
y = 0 a expansão é conhecida como Série de Maclaurin . Uma importante parte da expansão é a
forma polinomial:
Definição 2.4
Seja f (x) ∈ C n [a, b] e x0 ∈ [a, b]. O polinômio
n
X (x − x0 )k k
Tn [f, x0 ](x) = D f (x0 )
k=0
k!
é chamado de polinômio de Taylor de ordem n e centrado em x0 . ♣

Para provar alguns resultados nos próximos capı́tulos, vamos necessitar a troca da ordem da
integração dupla. Para isto, revemos o seguinte teorema,
Teorema 2.4 (Teorema de Fubini)
Sejam [a, b] e [c, d] dois intervalos compactos, f uma função Rieamann integravel e assumimos
quase Z b
g(y) = f (x, y)dx
a

existe para todo y fixo em [c, d]. Então g é Rieamann integravel em [c, d] e
Z Z d Z b
f (x, y)d(x, y) = f (x, y)dx dy. (2.1)
[a,b]×[c,d] c a

15
2.1 Diferenciação e Integração

Além disso, se Z d
h(y) = f (x, y)dy
c

existe para todo x fixo em [a, b], então

Z b Z d Z d Z b Z
f (x, y)dy dx = f (x, y)dx dy = f (x, y)d(x, y). (2.2)
a c c a [a,b]×[c,d]
♡

Abordamos a maioria dos resultados analı́ticos sobre diferenciação de ordem inteira e integração,
que vamos precisar em capı́tulos posteriores sobre Cálculo Fracionário ou reformular na configuração
fracionária. No entanto, para algumas provas nas próximas seções e capı́tulos vários teoremas de
ponto fixo serão necessários, então os enunciamos nesta seção, embora eles não são resultados no
campo da integração e diferenciação em sentido estrito. Iniciamos com Teorema do Ponto Fixo de
Weissinger. Algumas definições sobre espaços métricos:
Definição 2.5
Seja (E, d) espaço métrico e F ⊆ E.
Espaço Completo: F é completo se toda sequência de Cauchy de elementos de F converve em
F
Equicontinuo: F é equicontinuo se, para todo ϵ > 0, ∃δ > 0 tal que, ∀f ∈ F, x, y ∈ [a, b]
com |x − y| < δ temos |f (x) − f (y)| < ϵ
Uniformemente Limitado: F é uniformemente limitado se existe a constante C tal que
∥f ∥∞ ≤ C para toda f ∈ F .
Compacto: F é compacto se toda sequência em F tem uma subsequencia convergente emF ,
Relativamente Compacto: O conjunto F é dito ser relativamente compacto em E se o fecho
de F é um subconjunto compacto de E, ♣

Teorema 2.5 (Teorema do Ponto Fixo de Weissinger)

Assuma (U, d) um espaço métrico completo não vazio, e seja αj ≥ 0 para j ∈ N0 e tal que
X∞
αj converge. Além disso, seja a aplicação A : U → U satisfaz
j=0

d(Aj u, Aj v) ≤ αj d(u, v) (2.3)

para todo j ∈ N e todo u, v ∈ U Então A tem um único ponto fixo definido u∗ . Além disso,
para qualquer u0 ∈ U , a sequencia (Aj u0 )∞ ∗
j=1 converge para o ponto fixo u .
♡

Uma consequência imediata do Teorema de Weissinger é o Teorema do Ponto de Fixo de Banach.

Corolário 2.1 (Teorema do Ponto Fixo de Banach)
Assuma (U, d) um espaço métrico completo não vazio, e seja 0 ≤ α < 1 . Além disso, seja a
aplicação A : U → U satisfaz
d(Au, Av) ≤ αd(u, v) (2.4)

16
2.2 Equações Diferenciais

para todo u, v ∈ U . Então A tem um único ponto fixo definido u∗ . Além disso, para qualquer
u0 ∈ U , a sequencia (Aj u0 )∞ ∗
j=1 converge para o ponto fixo u .
♡

Além disso, precisaremos de um resultado ligeiramente diferente que afirme apenas a existência,
mas não a singularidade de um ponto fixo em seções posteriores deste texto. Uma hipótese maı́s fraca
sobre o operador em questão leva ao Teorema de Schauder.
Teorema 2.6 (Teorema do Ponto Fixo de Schauder)
Seja (E, d) um espaço métrico completo, e seja U um subespaço convexo e fechado de E. Além
disso, seja a aplicação A : U → U tal que o conjunto {Au : u ∈ U } é relativamente compacto
em E. Então A tem ao menos um ponto fixo. ♡

Outro resultado que será útil na nossa análise:

Teorema 2.7 (Arzelà-Ascoli)
Seja F ⊆ C[a, b] com a < b, e assuma a norma de Chebyshev esteja definida nestes conjuntos.
Então F é relativamente compacto em C[a, b] se e somente se F é equicontı́nuo e uniformemente
limitado. ♡

Estes são os resultados básicos do cálculo clássico que usaremos nos próximos capı́tulos. Na
seção seguinte iremos abordar alguns resultados de equações diferenciais ordinárias no contexto
analı́tico e numérico.

2.2 Equações Diferenciais

Nesta seção iremos fazer uma revisão da teoria das equações diferenciais ordinárias, no as-
pecto analı́tico. O objetivo é abordar resultados que serão importante na formulação para o Cálculo
Fracionário. Iniciamos com a definição de uma equação diferencial ordinária.
Definição 2.6 (Equação Diferencial Ordinária)
Seja n ∈ N e f : A ⊆ R2 → R. Então
Dn y(x) = f (x, y(x)) (2.5)
é chamada de equação diferencial ordinária de ordem n. Se condições iniciais da forma
Dk y(0) = bk (k = 0, 1, 2, . . . n − 1) (2.6)
são definidas, entendemos as equações (2.5) e (2.6) como um problema de valor inicial (P.V.I.).
♣

Um primeiro resultado fundamental na formulação do Cálculo Fracionário, é dado pelo seguinte

lema:

17
2.2 Equações Diferenciais

Lema 2.3
A função y(x) é solução do P.V.I (2.5) e (2.6) se e somente se y(x) é uma solução da equação
integral
n−1 k Z x
X x k 1
y(x) = D y(0) + (x − t)n−1 f (t, y(t))dt (2.7)
k=0
k! (n − 1)! 0

Demonstração A demonstração é obtida aplicando o clássico operador diferencial ou integral em

conexão com o Lema 2.1 para obter a condição necessária e a suficiente
As questões de existência e unicidade do P.V.I (2.5) e (2.6) são respondidas pelos teoremas que
seguem:
Teorema 2.8 (Teorema da Existência de Peano)
Seja c > 0 e G := {(x, y) ∈ R2 : 0 ≤ x ≤ c} e f : G → R uma função contı́nua. Entâo o
P.V.I. (2.5) e (2.6) tem ao menos uma solução na vizinhança U ⊆ [0, 0 + c]
♡

Demonstração Ver [32] pág. 45.

Teorema 2.9 (Teorema da Existência e Unicidade de Picard-Lindelöf)
Seja c > 0 e G := [0, 0 + c] × R e assuma que a função f : G → R é contı́nua e que satisfaz a
condição de Lipschitz com respeito a segunda variável, isto é, existe uma constante L > 0, tal
que, para todo (x, y1 ) e (x, y2 ) ∈ G nós temos
|f (x, y1 ) − f (x, y2 )| < L|y1 − y2 |.
Entâo o P.V.I. (2.5) e (2.6) possui uma única solução.
♡

Demonstração Ver [32] pág. 48.

Outro resultado importante é determinar de como o comportamento da função f influencia
no comportamento da solução da equação diferencial (2.5). Para vermos estes resultados, vamos
introduzir a definição:
Definição 2.7
Seja n ∈ N, G ⊆ Rn e f ∈ C(G). Então a função f é chamada de analı́tica em G, se para
quaisquer pontos (v1 , v2 , . . . , vn ) ∈ G existe uma série de potências satisfazendo
∞
X
f (x1 , x2 , . . . , xn ) = Cµ1 ,...,µn (x − v1 )µ1 (x − v2 )µ2 · · · (x − vn )µn
µ1 ,...,µn =0

a qual é absolutamente convergente na vizinhança de (v1 , v2 , . . . , vn ).

♣

O teorema que segue, estabelece as condições para que a solução da equação diferencial (2.5)
apresenta uma solução analı́tica.

18
2.3 Esquemas Numéricos

Teorema 2.10
Se a função f do P.V.I. (2.5) e (2.6) é uma função analı́tica na vizinhaça de
(0, Dy(0), . . . Dn y(0)) a solução de (2.5) é analı́tica na vizinhaça de 0.
♡

Demonstração O Teorema de Picard-Lindelöf garante a existência e unicidade de uma solução para

um P.V.I. em uma vizinhança de um ponto inicial quando a função f e a sua derivada parcial em
relação a y são contı́nuas. A função analı́tica f é, por definição, infinitamente diferenciável e pode
ser representada por uma série de Taylor. Se f é analı́tica, então as suas derivadas parciais também
são analı́ticas, o que garante a continuidade necessária para o Teorema de Picard-Lindelöf. A solução
do P.V.I., obtida através do Teorema de Picard-Lindelöf, é então também uma função contı́nua e
infinitamente diferenciável. Além disso, a solução pode ser expressa como uma série de Taylor, o que
a define como uma função analı́tica
A diferenciabilidade da solução está relacionada com a diferenciabilidade de f , como podemos
ver no teorema que segue.
Teorema 2.11
Seja k ∈ N, b > 0 e f ∈ C k ([0, b] × R). Então a solução do equação diferencial
Dy(x) = f (x, y(x)), y(0) = b0
é (k+1)-vezes diferenciável. ♡

Como exemplo vamos considerar o problema da Lei de resfriamento de Newton, representado de

valor inicial dado por 
 dT
= −k(T − Tm )
dt (2.8)
 T (0) = T
0

onde k é uma constante positiva que depende do material que constitui o corpo e Tm é a temperatura
ambiente considerada constante. A solução analı́tica se apresenta na forma
T (t) = Tm + (T0 − Tm )e−kt (2.9)
Na Figura 2.1 é apresentado as curvas integrais para valores de Tm = 27, k = 1 e T0 no intervalo
[0, 20].

2.3 Esquemas Numéricos

No contexto das soluções numéricas das equações diferenciais ordinárias de primeira ordem com
condição inicial, estamos interessados em solução em um intervalo fechado [0, X], para algum X > 0.
Em geral um método numérico fornece a solução num conjunto de pontos {x0 , x1 , . . . , xN ∈ [0, X]},
onde xk+1 = xk + h, x0 = 0 e h = X/N . Sendo ym a aproximação de y(xm ) e fm = f (xm , ym ) a
discretização da equação diferencial, podemos definir o método linear de multipassos:

19
2.3 Esquemas Numéricos

Figura 2.1: Solução do p.v.i. 2.8 para T0 ∈ [0, 20]

Definição 2.8
Definimos um método linear de multipassos para uma equação diferencial ordinária de primeira
ordem por
p p
X X
αk ym−k = h βk f (xm−k , ym−k ) (2.10)
k=−1 k=−1

onde αk , βk para k = −1, 0, 1, . . . , p denota constantes reais. ♣

A cada método linear de multipassos está associado ao primeiro e segundo polinômio carac-
terı́stico, sendo estes
p
X
ρ(ζ) = αk ζ p−k (2.11)
k=−1
p
X
σ(ζ) = βk ζ p−k (2.12)
k=−1
(2.13)
A estrutura do método será denotada pelo tipo dos polinômios caracterı́sticos. Dizemos que o método
será do tipo (ρ, σ). O Lema 2.2 nos garante que uma equação diferencial ordinária é equivalente a um
tipo especı́fico de equação integral de Volterra. Desta forma, vamos formular um método linear de
multipassos (ρ, σ) aplicado a uma equação integral.
Lema 2.4
Um método linear multipassos (ρ, σ) aplicado a equação integral
Z x
y(x) = f (t)dt,
0

20
2.3 Esquemas Numéricos

pode ser descrito como

m
X
(h If )(x) = h ωm−j f (jh), x = mh.
j=0
onde h I denota a discretização do operador integral I Os pesos ωm são dados pela série de
potência que gera a função
σ(1/ζ)
ω(ζ) = (2.14)
ρ(1/ζ)
♡

Demonstração Este Lema é um caso particular do Lema 2.1 em [27]

Um método linear multipassos (ρ, σ), com p-passos necessita que seja conhecida a solução nos
pontos yi , i = 0, 1, . . . , p − 1, ou a solução não é obtida de forma única. Em [14] as questões de
existência, unicidade e convergência da aproximação da equação diferença para a solução da equação
diferencial são abordadas.
Definição 2.9
Seja f (x, y) definida em todo x ∈ [0, X] tal que o P.V.I.
Dy(x) = f (x, y(x)); y(0) = b0
tem solução para todo b0 . Um método linear de multipassos (ρ, σ) é dito ser convergente se
lim ym = y(xm )
h→0
para todo xm ∈ [0, X] e toda solução {ym } da equação diferença (2.10) temdo valores iniciais
y0 , . . . yp−1 satisfazendo
lim yi = b0 , i = 0, 1, . . . , p − 1
h→0
♣

Estabilidade e consistência de um método linear de multipassos equivale na convergência do

método. Este resultado é conhecido como Teorema de Equivalência de Lax [25]. (completar com
definição de estabilidade consistência e convergência, mesmo que seja no apendice)
Definição 2.10 (Estabilidade)
Um método linear de multipassos é estável, se e somente se, os correspondentes pesos ωm são
limitados. ♣

Definição 2.11 (Consistência)

Um método linear de multipassos é consistente de ordem p, se satisfaz
hω(e−h ) = 1 + O(hp ).
♣

Para finalizar este capı́tulo, apresentados dois exemplos de métodos linear de multipassos
Vamos considerar dois métodos que pertencem a famı́lia de métodos do tipo Adams, um explicito,
Adams-Bashforth e outro implı́cito, Adams-Moulton. Em ambos o primeiro polinômio caracterı́stico
é dado por
ρ(ζ) = ζ p+1 − ζ p

21
2.3 Esquemas Numéricos

a diferença entre eles está no segundo polinômio caracterı́stico, fornecendo duas diferentes equações
diferença:
p
X
ym+1 = ym + h βk f (tm−k , ym−k ) (Adams-Bashforth)
k=0
p
X
ym+1 = ym + h βk f (tm−k , ym−k ) (Adams-Moulton)
k=−1

onde os coeficientes βk são obtidos pela interpolação polinomial sobre os pontos tm , tm−1 , . . . tm−p ,
no caso do método de Adams-Bashforth e tm+1 , tm , . . . tm−p , no caso do método de Adams-Moulton.
(p)
As Tabelas 2.1 e 2.2 apresentam os coeficientes βk , onde p determina o número de passos do
método.
Tabela 2.1: Coeficientes de Adams-Bashforth da segunda caracterı́stica polinomial
k 0 1 2 3 4 5
(0)
βk 1
(1)
2 βk 3 -1
(2)
12 βk 23 -16 5
(3)
24 βk 55 -59 37 -9
(4)
720 βk 1901 -2774 2616 -1274 251
(5)
1440 βk 4227 -7673 9482 –6798 2627 -425

Tabela 2.2: Coeficientes de Adams-Moulton da segunda caracterı́stica polinomial

k -1 0 1 2 3 4
(0)
βk 1
(1)
2 βk 1 1
(2)
12 βk 5 8 -1
(3)
24 βk 9 19 -5 1
(4)
720 βk 251 646 -264 106 -19
(5)
1440 βk 475 1427 -798 482 -173 27

Outra classe de esquemas numérico pode ser obtido na aproximação da derivada da função por
diferenças finitas. Consideremos Esquemas de Diferenças Finitas de passo múltiplo de ordem p, cuja
a fórmula geral é dada por:
p
X
αk ym−k = hf (xm+1 , ym+1 )
k=−1

onde ym−k representa a aproximação da solução na etapa m − k, h é o tamanho do passo, αk são os

coeficientes do método e f (xm+1 , ym+1 ) é a função que define a equação diferencial avaliada na etapa
m + 1.
Os coeficientes αk , de ordem p, podem ser descritos como os coeficientes de uma série de
Maclaurin de uma função geradora correspondente. A função geradora para os coeficientes αk é dada

22
2.4 Exemplos Numéricos

por:
p p
X
k
X 1
α(z) = αk z = (1 − z)k
k=0 k=1
k
(p)
A Tabela 2.3 mostra os coeficientes αk . Os coeficientes são calculados com base na função
geradora mencionada acima.

Tabela 2.3: Coeficientes da Fórmula de Diferenças Finitas da primeira caracterı́stica polinomial

k -1 0 1 2 3 4 5
(0)
αk 1
(1)
αk 1 -1
(2)
2 αk 3 -4 1
(3)
6 αk 11 -18 9 -2
(4)
12 αk 25 -48 48 16 3
(5)
60 αk 137 -300 300 -200 75 -12
(6)
60 αk 147 -360 450 -400 225 -72 10

2.4 Exemplos Numéricos

Iremos testar os esquemas numéricos aplicados ao problema da Lei de resfriamento de Newton
2.8, com T0 = 0, 5, 10, Tm = 27 e k = 1.
Iniciamos com os Esquemas Numéricos de Adams-Bashforth de passos p = 1, 2, 3, 4, onde os
valores iniciais y1 , . . . yp−1 são obtidos pelo Esquema Numérico de Euler Melhorado. (Ver Apêndice
A).
Testamos o esquema de passo um, usando uma malha de 1024 pontos para aproximar a solução
do problema de valor inicial 2.8, com T0 = 0, 10, 15, 20. Neste caso, necessitamos calcular uma
aproximação y1 = y(t1 ) e neste caso foi utilizado o Esquema de Euler Melhorado. As aproximações
são apresentadas na Figura 2.2.
Na Figura 2.3 temos o erro absoluto E = ∥ym −y(xm )∥ com uma malha de 1024 pontos, T0 = 10
e p = 1, 2, 3, 4. O que podemos observar é que existe uma diferença significativa entre o esquema
de passo um com os demais. Em todos os casos o erro tende a zero, devido ao comportamento do
problema, que tende a temperatura ambiente, que neste exemplo temos Tm = 27.
Também testamos um esquema do tipo Preditor-Corretor, onde o Esquema de Euler Melhorado
atua como preditor e o Esquema de Adams-Moulton como corretor. A implementação está descrita
em Apêndice B. Testamos a precisão dos Esquemas de Adams-Moulton para aproximar a solução do
problema de valor inicial 2.8, com T0 = 0, Tm = 27, k = 1.
Na Figura 2.4 temos o erro absoluto E = ∥ym − y(xm )∥ com uma malha de 1024 pontos, T0 = 0
e p = 1, 2, 3, 4. O que podemos observar é que existe uma diferença significativa entre o esquema de
passo um com os demais, como também observado nos esquemas explı́citos. Além disso, a ordem do
erro é menor que o erro obtido com os esquemas explicitos.

23
2.4 Exemplos Numéricos

Figura 2.2: Aproximação Adams-Bashforth com N = 1024, no intervalo de tempo [0, 10].

Figura 2.3: Aproximação com N = 1024, no intervalo de tempo [0, 10], T0 = 0 e p = 1, 2, 3, 4

24
2.4 Exemplos Numéricos

Figura 2.4: Preditor-Corretor com N = 1024, no intervalo de tempo [0, 10], T0 = 0 e p = 1, 2, 3, 4

25
Capı́tulo 3 Transformadas Integrais e Funções
Especiais
Neste capı́tulo iremos fazer uma revisão de algumas transformadas integrais e funções especiais,
que serão essenciais na definição do Cálculo Fracionário. Iniciamos o capı́tulo, abordando a Trans-
formada de Fourier. Na solução de equações diferenciais lineares, a Transformada de Fourier pode
ser usada para transformar a equação para o domı́nio de frequência, onde pode ser resolvida com
métodos algebraicos. Depois, a solução pode ser transformada de volta ao domı́nio do tempo com
a utilização da Transformada de Fourier Inversa. Também iremos ver a Transformada de Laplace,
que no contexto de solução de equações diferenciais lineares, a Transformada de Laplace pode ser
usada para transformar a equação para o domı́nio de Laplace, onde pode ser resolvida com métodos
algébricos. Depois, a solução pode ser transformada de volta ao domı́nio do tempo com a utilização
da Transformada de Laplace inversa.
As funções Gama de Euler e a função de Mittag-Leffler, generalizam situações que se são
importante no desenvolvimentos do Cálculo Fracionário. Por essa razão, faremos uma revisão das
definições e propriedades que serão essenciais para o próximo capı́tulo.

3.1 Transformadas Integrais

Iniciamos com a definição de uma transformada integral.
Definição 3.1
Seja a função φ(x) pertencente a um espaço de funções. Então temos uma transformada integral
dada por Z +∞
(Kφ)(x) = k(x, t)φ(t)dt = g(x), (3.1)
−∞

onde k(x, t) é alguma função, chamada de kernel da transformação, e g é a transformação da

função φ.
♣

Iniciamos com a Transformada de Fourier unidimensional, onde a função kernel é da forma

k(x, t) = eixt , com x ∈ R.
Definição 3.2 (Transformada de Fourier)
Seja a função φ(x) definida para x ∈ R, a transformada integral
Z +∞
(Fφ)(x) = eixt φ(t)dt = φ̂(x) (3.2)
−∞
é chamada de Transforma da Fourier. cuja a transformada inversa é dada por
Z +∞
−1 1 1
(F g)(x) = e−ixt g(t)dt = ĝ(−x) (3.3)
2π −∞ 2π
♣
3.1 Transformadas Integrais

Sendo as Transformadas 3.2 e 3.3 inversas uma da outra, ou seja, (FF −1 φ)(x) = (F −1 Fφ)(x) =
φ(x), segue que
(FFφ)(x) = φ(−x) (3.4)

Outras propriedades bem conhecidas da Transformada de Fourier são

F[Dn φ](x) = (−ix)n (Fφ)(x), n ∈ N (3.5)
e
Dn [Fφ](x) = (ix)n (Fφ)(x), n ∈ N (3.6)

onde Dn denota o clássico operador diferencial de ordem n.

No caso n-dimensional a Transformada de Fourier é definida da seguinte forma:
Definição 3.3 (Transformada de Fourier em Rn )
Seja a função φ(x) definida para x ∈ Rn , a transformada integral
Z
(Fφ)(x) = eixt φ(t)dt = φ̂(x) (3.7)
Rn
cuja a transformada inversa é dada por
Z
−1 1 1
(F g)(x) = n
e−ixt g(t)dt = ĝ(−x) (3.8)
(2π) Rn (2π)n
♣

As propriedades apresentadas no caso unidimensional também são válidas para o caso n-

dimensional.
No caso da Transformada de Laplace, a função kernel é da forma
(
e−xt se t > 0,
k(x, t) =
0 se t ≤ 0.
Com isto temos a seguinte definição
Definição 3.4 (Transformada de Laplace)
Seja a função φ(t) definida para t ∈ (0, ∞), a transformada integral
Z +∞
(Lφ)(s) = e−st φ(t)dt = φ̃(s) (3.9)
0
é chamada de Transforma da Laplace. cuja a transformada inversa é dada por
Z γ+bi
−1 1
(L φ)(s) = lim e−ixt φ(p)dp (3.10)
2πi b→∞ γ−bi
onde a integração é feito ao longo da linha vertical Re(p) = γ no plano complexo em que γ é
maior do que a parte real de todas as singularidades de φ(p).
♣

Algumas propriedades da Transformada de Laplace são análogas as apresentadas para a Trans-

formada de Fourier.
As transformadas 3.9 e 3.10 são inversas uma da outra, ou seja, (LL−1 φ)(x) = (L−1 Lφ)(x) =
φ(x), quando φ é uma função “suficientemente boa”.

27
3.2 Função Gama de Euler

As propriedades relacionadas com o operador diferencial são dadas por:

n−1
X
n n
L[D φ](s) = s (Lφ)(s) − sn−j−1 (Dj φ)(0), n ∈ N (3.11)
j=0
e
Dn [Lφ](s) = (−1)n (Ltn φ(t))(s), n ∈ N (3.12)

onde Dn denota o clássico operador diferencial de ordem n.

A ideia da convolução é permitir uma análise de sistemas lineares e invariantes no tempo. Isto é,
sistemas que têm uma resposta proporcional à entrada, sem se alterar com o passar do tempo.
Nas transformadas integrais, por exemplo, a convolução de duas funções no tempo é equivalente
ao produto de suas transformadas. Definimos a convolução de duas funções por
Definição 3.5 (Convolução)
Sejam as funções f, g ∈ L1 (R). Definimos a convolução de f com g por:
Z ∞
(f ∗ g)(x) := f (x − t)g(t)dt.
−∞
♣

A convolução de duas funções satisfazem as seguintes propriedades:

Teorema 3.1
Sejam as funções f, g, h ∈ L1 (R) e a ∈ C. Então as seguintes propriedades são válidas:
1. Comutativa: f ∗ g = g ∗ f ,
2. Associativa: f ∗ (g ∗ h) = (f ∗ g) ∗ h,
3. Distributiva: f ∗ (g + h) = f ∗ g + f ∗ h,
4. Associativa com Escalar: a(f ∗ g) = (af ) ∗ g = f ∗ (ag)
♡

A transformada integral da convolução de duas funções no tempo equivalente ao produto de suas

transformadas. No caso das Transformadas de Fourier e Laplace, temos.
√
F(f ∗ g)(x) = 2π(Ff )(x) · (Fg)(x) e L(f ∗ g)(x) = (Lf )(x) · (Lg)(x)

3.2 Função Gama de Euler

Iniciamos o estudo de funções especiais, pela Função Gama de Euler. Esta função pode ser
encontrada em vários tópicos, em especial na definição de outras funções especiais de forma direta ou
por integração de contorno. Iniciamos com a definição de Constante de Euler-Mascheroni.
Definição 3.6 (Constante de Euler-Mascheroni)
A Constante de Euler γ é dada por
n
X 1
γ = lim − ln(n) = 0.577215664... (3.13)
n→∞
k=1
k
♣

Existem várias formas de se definir a Função Gama de Euler. Adotamos aquela que será mais

28
3.2 Função Gama de Euler

Figura 3.1: A Função Gama (vermelha) e sua reciproca (azul).

útil no contexto do Cálculo Fracionário.

Definição 3.7 (Função Gama de Euler)
Para z ∈ C\Z− , definimos a função Gama de Euler por
 Z +∞
tz−1 e−t dt se Re(z) > 0




0
Γ (z) = (3.14)



 Γ (z + 1)/z se Re(z) ≤ 0,
♣

A Função Gama de Euler é uma função meromorfa, isto é, uma função complexa diferenciável
em todo o plano complexo exceto em um número de pontos, chamados de polos. Neste caso os polos
são os inteiros negativos, Z− .
Na Figura 3.1 temos a representação da Função Gama de Euler e sua reciproca. Veremos algumas
propriedades que a Função Gama satisfaz.
Teorema 3.2
A Função Gama satisfaz as seguintes propriedades:
1. Para Re(z) > 0 segue que
Z 1 z−1
1
Γ (z) = ln dt
0 t
2. Para z ∈ C\Z−
Γ (z + 1) = zΓ (z)

3. Para n ∈ N
Γ (n) = (n − 1)!

4. Para z ∈ C\Z−
Γ (1 − z) = −zΓ (−z)

29
3.2 Função Gama de Euler

5. Para Re(z) > 0 temos a representação por limite

n!nz
Γ (z) = lim
n→∞ (z + 1)(z + 2) · · · (z + n)

6. Para z ∈ C\Z− segue a definição de Weierstrass

∞
1 γz
Y z −z/n
= ze 1+ e
Γ (z) n=1
n
onde γ é a constante de Euler (3.13).
7. A Função Gama é analı́tica para todo z ∈ C\Z−
8. A Função Gama não se anula.
9. (Reflexão) Para todo z ∈ C não inteiro
π π
Γ (z)Γ (1 − z) = e Γ (z)Γ (−z) = −
sen(πz) zsen(πz)
10. Para n ∈ N, e n ı́mpar.
n (n − 2)!!√π
Γ =
2 2(n−1)/2
onde n!! é o duplo fatorial definido por

 n · (n − 2) · · · 5 · 3 · 1 n > 0 e ı́mpar

n!! = n · (n − 2) · · · 6 · 4 · 2 n > 0 e par

1 n = 0, −1

♡

Demonstração
1. Fazendo a troca de variáveis u = −log(1/t) temos que du = 1/tdt e t = −e−u Quanto t = 1 ⇒
u = 0 e quando t → 0 ⇒ u → ∞. Logo obtemos que
Z 1 z−1 Z 0 Z +∞
1 −u
ln dt = z−1
u (−e )du = uz−1 e−u du
0 t +∞ 0
2. Aplicando a integração por partes na integral de (3.14 ), onde u = tz−1 ⇒ du = (z − 1)tz−1 e
dv = e−t ⇒ v = −e−t temos Z ∞
Γ (z) = tz−1 e−t dt
0
Z ∞
∞
= −tz−1 e−t 0 + (z − 1)tz−1 e−t dt
0
Z ∞
z−1 −t
= (z − 1) t e dt
0
= (z − 1)Γ (z − 1)
3. Como Γ (1) = 1 e do resultado da propriedade 2. segue que
Γ (n) = (n − 1)Γ (n − 1) = (n − 1)(n − 2)Γ (n − 2) = . . . = (n − 1)(n − 2) · 2Γ (1) = (n − 1)!
4. Da propriedade 2. temos que se Re(−z) > 0 segue
Γ (1 − z) = Γ (1 + (−z)) = −zΓ (−z)

30
3.2 Função Gama de Euler

G Quando Re(−z) ≤ 0 da definição (3.14) segue que

Γ (−z + 1)
Γ (−z) = ⇒ Γ (1 − z) = −zΓ (−z)
−z
5. Para a demostração introduzimos a função auxiliar
Z n n
t
Γn (z) = 1− tz−1 dt
0 n
Primeiramente faremos uma mudança de variável da forma s = t/n ⇒ ds = 1/ndt. Quando
t = 0 ⇒ s = 0 e quando t = n ⇒ s = 1. Logo
Z 1 Z 1
n z−1 z
Γn (z) = (1 − s) (ns) nds = n (1 − s)n sz−1 ds
0 0
Agora aplicamos a integração por partes n − vezes, obtendo
n 1 nz n 1
Z Z
z n n 1 ( z
Γn (z) = n (1 − s) s |0 + (1 − s) n − 1)s ds = (1 − s)( n − 1)sz ds
z 0 z 0
Z 1
z
nz n(n − 1) 1

n−1
Z
n n−1 n+1 1 ( z+1
= n (1 − s) s 0
+ (1 − s) n − 2)s ds = (1 − s)( n − 1)
z z+1 0 z(z + 1) 0
Z 1
nz n!
= sz+n−1 ds
z(z + 1) · · · (z + n − 1) 0
nz n!
=
z(z + 1) · · · (z + n)

Observe que Γ (z) = limn→∞ Γn (z) o que conclui da demostração da propriedade.

6. A função auxiliar também será usada para esta demonstração. Observe que
nz n!
Γn (z) =
z(z + 1) · · · (z + n)
nz
= .
z(1 + z/1) · · · (1 + z/n)
Usando a relação
nz = ez ln(n) = ez(ln(z)−1−1/2−...−1/n) ez+z/2+...+z/n

segue que
nz
Γn (z) =
z(1 + z/1) · · · (1 + z/n)
ez(ln(z)−1−1/2−...−1/n) ez+z/2+...+z/n
=
z(1 + z/1) · · · (1 + z/n)
1 ez ez/2
= ez(ln(z)−1−1/2−...−1/n) . . . ez/n (1 + z/n).
z (1 + z) (1 + z/2)

31
3.2 Função Gama de Euler

Com a representação por limite da Função Gama da propriedade 5. segue que

1 1
= lim
Γ (z) n→∞ Γn (z)

= ezγ z lim e−z (1 + z)e−z/2 (1 + z/2) . . . e−z/n (1 + z/n)

n→∞
∞
γz
Y z −z/n
= ze 1+ e
n=1
n
7. Segundo a definição de Weierstrass, uma função é analı́tica para todo z finito e seus únicos zeros
são simples em z = 0 e em cada inteiro negativo. Assim da propriedade 6. segue que a Função
Gama é analı́tica.
8. Da propriedade 6. podemos ver que 1/Γ (z) não possui polos. Logo Γ (z) não se anula.
9. Da propriedade 6. segue que
∞ h
1 1 2 γz −γz
Y z −z/n z z/n i
= z e e 1+ e 1− e
Γ (z) Γ (−z) n=1
n n
∞
z2
Y
2
= −z 1− 2
n=1
n
Da propriedade 4. substituirmos Γ (−z) por Γ (1 − z)/(−z), obtendo
∞
z2

1 1 Y
=z 1− 2
Γ (z) Γ (1 − z) n=1
n
e comparando com a série
∞
z2
Y
sen(πz) = πz 1− 2
n=1
n
chegamos ao resultado
π
Γ (z)Γ (1 − z) =
sen(πz)
A segunda parte da propriedade segue de forma análoga, se fazer a substituição do resultado da
propriedade 4.
10. Sendo n ∈ N e n ı́mpar e da propriedade 2. seque que

n n−2+2 n−2 n−2 n−2
Γ =Γ =Γ +1 = Γ
2 2 2 2 2
Repetindo o processo seque
n (n − 2) (n − 4) (n − 6) √
1 1 (n − 2)!! π
Γ = ... Γ =
2 2 2 2 2 2 2(n−1)/2

Muitas outras propriedades da Função Gama de Euler podem ser encontradas em [3]. As
apresentadas no Teorema (3.2) serão suficiente para os próximos capı́tulos. A propriedade 3. mostra
que a Função Gama é uma generalização do fatorial, que nos leva a uma generalização dos coeficientes
binomiais para α ∈ R

32
3.3 A Função Beta

Definição 3.8
Para α ∈ R e k ∈ N∗ definimos
(−1)k−1 Γ (k − α)

α α(α − 1) · · · (α − k + 1)
= =
k Γ (1 − α)Γ (k + 1) k!
♣

3.3 A Função Beta

A Função Beta é uma função especial de duas variáveis conectada diretamente com a Função
Gama, definida da seguinte forma:
Definição 3.9 (Função Beta)
A Função Beta B(z, w), com z, w ∈ C é definida por:
Γ (z)Γ (w)
B(z, w) = (3.15)
Γ (z + w)
♣

Vamos estabelecer algumas propriedades da função, as quais serão necessárias para estabelecer
resultados nos próximos capı́tulos, especialmente a integral de Beta, que será usada nos exemplos no
Cálculo Fracionário. .
Teorema 3.3
A Função Beta satisfaz as seguintes propriedades:
1. Para Re(z), Re(w) > 0 (3.15) é equivalente a
Z 1
B(z, w) = tz−1 (1 − t)w−1 dt
Z0 +∞
tz−1
= dt
0 (1 + t)z+w
Z π/2
= 2 (sen t)2z−1 (cos t)2w−1 dt
0
2. A função B(z + 1, w + 1) é solução da Beta Integral:
Z 1
tz (1 − t)w dt = B(z + 1, w + 1)
0
3. B(z, w) = B(w, z)
4. B(z, w) = B(z + 1, w) + B(z, w + 1)
w w
5. B(z, w + 1) = B(z + 1, w) = B(z, w)
z z+w ♡

Demonstração
1. Primeiramente, vamos verificar as equivalências entre as integrais. Aplicando a mudança de
variável t = x/(x + 1) segue que dt = dx/(x + 1)2 , logo
Z 1 Z ∞ z−1 w−1
z−1 w−1 x x
t (1 − t) dt = 1− (x + 1)2 dx
0 x+1 x+1
Z0 ∞ z−1
x
= dx
0 (1 + x)z+w

33
3.3 A Função Beta

Reescrevendo a integral com a variável t, obtemos a primeira igualdade. Para a igualdade com
a última integral fazemos a mudança de variável t = sen2 (φ) e temos dt = 2sen(φ) cos(φ)dφ.
Z 1 Z π/2
z−1 w−1
t (1 − t) dt = (sen2 φ)z−1 (1 − sen2 φ)w−1 2senφ cos φdφ
0 0
Z π/2
= 2 (senφ)2z−1 (cos φ)2w−1 dφ
0
Novamente, reescrevendo na variável t obtemos o resultado. Agora vamos mostrar que essas
integrais são equivalentes a Função Beta. Para isto iniciamos com a relação.
Z ∞ Z ∞
−t z−1
Γ (z)Γ (w) = e t dt e−s sw−1 ds
0 0
Fazendo as mudanças de variáveis t = x2 ⇒ dt = 2xdx e s = y 2 ⇒ ds = 2ydy segue
Z ∞ Z ∞
−x2 2z−2 2
Γ (z)Γ (w) = 4 e x xdx e−y y 2w−2 ydy
Z0 ∞ Z ∞ 0
2 2
= 4 e−x −y x2z−1 y 2w−1 dxdy
0 0

Usando coordenadas polares x = rcosθ e y = rsenθ reescrevemos a equação acima da forma

Z ∞ Z π/2
2
Γ (z)Γ (w) = 4 e−r r2z+2w−2 (cos θ)2z−1 (senθ)2w−1 rdθdr
0 0
Z ∞ Z Z π/2
−r2 2z+2w−1
= 2 e r dr × 2 (cos θ)2z−1 (senθ)2w−1 dθ
0 0

√ √
Com r = t ⇒ dr = dt/(2 t) e θ = π/2 − φ finalmente chegamos que
Z ∞ Z Z π/2
−t z+w−1
Γ (z)Γ (w) = e t dt × 2 (senφ)2z−1 (cos φ)2w−1 dφ
0 0
Z π/2
= Γ (z + w) × 2 (senφ)2z−1 (cos φ)2w−1 dφ
0

Ao compararmos com a definição da Função Beta B(z, w), concluı́mos a propriedade.

2. Basta considerar a equivalência com a primeira integral da propriedade 1.
3. Da equivalência com a primeira integral da propriedade 1.
Z 1
B(z, w) = tz−1 (1 − t)w−1 dt
Z0 0
= (1 − u)z−1 (u)w−1 − du
1
Z 1
= 2 uw−1 (1 − u)z−1 du
0
= B(w, z)
onde foi aplicado a mudança de variável u = 1 − t.

34
3.4 Função de Mittag-Leffler

4. Da equivalência com a segunda integral da propriedade 1.

Z +∞
tz−1
B(z, w) = dt
0 (1 + t)z+w
Z +∞
t+1 tz−1
= dt
0 t + 1 (1 + t)z+w
Z +∞ Z +∞
t tz−1 1 tz−1
= dt + dt
0 t + 1 (1 + t)z+w 0 t + 1 (1 + t)z+w
Z +∞ Z +∞
tz tz−1
= dt + dt
0 (1 + t)z+w+1 0 (1 + t)z+w+1
= B(z + 1, w) + B(z, w + 1)
5. Para a primeira igualdade aplicamos a integração por partes na primeira equivalência da propriedade
1. com u = (1 − t)w ⇒ du = −w(1 − t)w−1 dt e dv = tz−1 ⇒ v = tz /z sendo.
Z 1
B(z, w + 1) = tz−1 (1 − t)w dt
0
z 1
w 1 z
Z
wt
= (1 − t) + t (1 − t)w−1 dt
z 0 z 0
w
= B(z + 1, w)
z

Para provar a segunda igualdade vamos usar o resultado acima e da propriedade 4. Do resultado da
primeira igualdade temos
z
B(z + 1, w) = B(z, w + 1)
w
Substituindo no resultado da propriedade 4. segue
B(z, w) = B(z + 1, w) + B(z, w + 1)
z
= B(z, w + 1) + B(z, w + 1)
wz
= + 1 B(z, w + 1)
w
z+w
= B(z, w + 1)
w

3.4 Função de Mittag-Leffler

Nesta seção iremos abordar a Função de Mittag-Leffler, que está conectada com a Função de
Euler e será útil na análise de equações diferenciais ordinárias.

35
3.4 Função de Mittag-Leffler

Definição 3.10 (Função de Mittag-Leffler)

Para z ∈ C definimos a Função de Mittag-Leffler Eα (z) por
+∞
X zk
Eα (z) = ,α>0 (3.16)
k=0
Γ (αk + 1)
e sua forma generalizada Eα,β (z) por
+∞
X zk
Eα,β (z) = , α, β > 0 (3.17)
k=0
Γ (αk + β)
♣

No próximo teorema destacamos algumas das propriedades da Função de Mittag-Leffler.

Teorema 3.4
A Função de Mittag-Leffler possui as seguintes propriedades:
1. Para |z| < 1 a função generalizada de Mittag-Leffler satisfaz.
Z +∞
1
e−t tβ−1 Eα,β (tα z)dt =
0 z−1
2. Para |z| < 1 a Transformada de Laplace da Função Mittag-Leffler tβ−1 Eα,β (ztα ) é dada
por
sα−β
L[tβ−1 Eα,β (atα )](s) = α
s −a
3. Para valores especiais de α temos
1
(a) E0 (z) = 1−z (b) E1 (z) = ez
(c) E2 (z 2 ) = cosh(z) (d) E2 (−z 2 ) = cos(z)
♡

Demonstração
1. Para esta demonstração iremos usar a Transformada de Laplace de tk eqt , sendo
k!
L(tk eqt ) = , Re(p) > |q|.
(q − p)k+1
Calculando a Transformada de Laplace da função Mittag-Leffler segue
Z +∞ Z +∞ +∞
−t β−1
X (ztα )k
α
e t Eα,β (zt )dt = e−t tβ−1 dt
0 0 k=0
Γ (αk + β)
∞ k Z +∞
X z
= e−t tαk+β−1 dt
k=0
Γ (αk + β) 0
+∞
X zk
= Γ (αk + β)
k=0
Γ (αk + β)
+∞
X 1
= zk = , |z| < 1
k=0
1−z

36
3.4 Função de Mittag-Leffler

Figura 3.2: Função Mittag-Leffler, com α = 1 (sólida), α = 2 (tracejada), α = 3 (traço-ponto) e

α = 4 (pontilhada)

2. Da definição da Transformada de Laplace, seque que

Z +∞
β−1 α
L[t Eα,β (at )](s) = e−st tβ−1 Eα,β (atα )dt
0
+∞ ∞
(atα )k
Z X
−st β−1
= e t dt
0 k=0
Γ (αk + β)
+∞ ∞
ak
X Z
= e−st tαk+β−1 dt
k=0
Γ (αk + β) 0
∞
X ak Γ (αk + β)
=
k=0
Γ (αk + β) sαk+β
∞
1 X a k
=
sβ k=0
sα
1 1
=
s 1 − a/sα
β

sα−β
= α
s −a
3. Para cada caso é só observar que a série de potências representa a função em questão.

Na Figura 3.2 apresentamos o gráfico de Função Mittag-Leffler, geradas pelo Octave e truncando
a soma em n = 20.

37
Capı́tulo 4 Cálculo Fracionário
Neste capı́tulo iremos abordar as formulações clássicas para a derivada e integral não inteira. Os
resultados desse capı́tulo pode ser encontrado nos livros [Oldhann74], [30], [1] e [45]. Os aspectos
da integração e derivação e equações diferenciais serão visto separadamente.
No Capı́tulo 1, fizemos uma breve revisão histórica e podemos concluir que não existe uma forma
única de transferir os conceitos do cálculo inteiro para o não inteiro. Uma vez que uma derivada
fracionária é uma generalização da derivada comum e perderá muitas de suas propriedades básicas.
Como exemplo, não temos uma interpretação geométrica ou fı́sica clara da derivada fracionária. A
lei do ı́ndice só é válido quando se trabalha em espaços funcionais especı́ficos, a derivada do produto
de duas funções é difı́cil de calcular e a regra da cadeia não pode ser aplicado diretamente. Logo é
natural perguntar, quais propriedades de derivadas fracionárias as tornam tão adequadas para modelar
certos sistemas complexos.
No contexto em que a integração e diferenciação são operadores funcionais, iremos dividir este
capı́tulo seções dedicadas a cada um dos operadores não inteiros e suas propriedades. Iniciamos com
o operador de Riemann-Liouville.

4.1 Operador de Riemann-Liouville

No Lema 2.1 apresentamos o operador integral no conceito clássico. Trocando n ∈ N por α ∈ R
e usando a Função Gamma (3.14) no lugar do fatorial, obtemos a seguinte versão do operador integral
não inteiro:
Definição 4.1 (Riemann-Liouville)
Seja α ∈ R+ . O operador Iaα definido em L1 [a, b] por
Z x
α 1
Ia f (x) = (x − t)α−1 f (t)dt, (4.1)
Γ (α) a
para a ≤ x ≤ b, é chamado de operador integral não inteiro de Riemann-Liouville. Quando
α = 0 definimos o operador identidade Ia0 := I.
♣

Nos casos em que α ∈ N, o operador de Riemann-Liouville coincide com o operador clássico

definido no Lema 2.1, com a exceção que o domı́nio é estendido de Funções Riemann integráveis para
funções Lebesgue integráveis. Pelo fato da integral ter o produto de uma função integrável e uma
função contı́nua, a existência da integral é assegurada para o caso de α > 1. A existência para o caso
de α ∈ (0, 1) é assegurada pelo seguinte teorema.
Teorema 4.1
Seja f ∈ L1 [a, b] e α > 0. Então a integral Iaα existe em quase todo x ∈ [a, b] e a função
Iaα f ∈ L1 [a, b].
♡
4.1 Operador de Riemann-Liouville

Demonstração Usamos as funções de ’ligar’ e ’desligar’ para reescrever a integral

Z x Z +∞
α−1
(x − t) f (t)dt = ϕ1 (x − t)ϕ2 (t)dt
a −∞
onde, (
uα−1 para 0 < u ≤ b − a
ϕ1 (u) =
0 caso contrário
e (
f (u) para 0 ≤ u ≤ b
ϕ2 (u) =
0 caso contrário
Pela construção, ϕj ∈ L1 (R) com j = 1, 2, e pelo resultado em integração de Lebesgue para a
convolução das funções, ϕ1 ∗ ϕ2 , o resultado do teorema é satisfeito
A seguir iremos mostrar algumas propriedades, satisfeitas pelo operador clássico, e que também
são satisfeitas pelo operador não inteiro de Riemann-Liouville.
Teorema 4.2
Os operadores {Iaα : L1 [a, b] → L1 [a, b]; α ≥ 0} formam um semigrupo comutativo com
respeito a concatenação. O operador Ia0 é o elemento neutro deste semigrupo. Isto implica
diretamente que:
Iaα Iaβ = Iaα+β = Iaβ+α (4.2)

se f ∈ L1 [a, b] e α, β ∈ R+ .
♡

Demonstração Para valores de α, β ∈ N este resultado é bem conhecido. O fato de que Ia0 ser
o elemento neutro do semigrupo decorre da Definição 4.1. Para a demostração de (4.2) usamos a
definição do operador não inteiro.
Z x Z t
α β 1 α−1
Ia Ia f (x) = (x − t) (t − τ )β−1 f (τ )dτ dt
Γ (α)Γ (β) a a
Com a existência das integrais garantida pelo Teorema 4.1 e pelo Teorema de Fubini 2.1 podemos
trocar a ordem de integração, obtendo
Z xZ x
α β 1
Ia Ia f (x) = (x − t)α−1 (t − τ )β−1 f (τ )dtdτ
Γ (α)Γ (β) a τ
Z x Z x
1
= f (τ ) (x − t)α−1 (t − τ )β−1 dtdτ
Γ (α)Γ (β) a τ

Fazendo a substituição da forma t = τ + s(x − τ ) segue que

Z x Z 1
α β 1
Ia Ia f (x) = f (τ ) [(x − τ )(1 − s)]α−1 × [s(x − τ )]β−1 (x − τ )dsdτ
Γ (α)Γ (β) a 0
Z x Z 1
1 α+β−1
= f (τ )(x − τ ) (1 − s)α−1 sβ−1 dsdτ
Γ (α)Γ (β) a 0

A segunda integral é a Função Beta (ver Definição 3.9 )

Γ (α)Γ (β)
B(α, β) = .
Γ (α + β)

39
4.1 Operador de Riemann-Liouville

Com isto obtemos Z x

1
Iaα Iaβ f (x) = f (τ )(x − τ )α+β−1 dτ = Iaα+β
Γ (α + β) a
válido em quase todo [a, b]
Os próximos teoremas tratam de propriedades de mapeamento dos operadores integrais não
inteiros, que serão importantes para o enunciações de outros teoremas.
Teorema 4.3
Seja f ∈ Hµ [a, b] para algum µ ∈ [0, 1] e seja 0 < α < 1. Então
f (a)
Iaα f (x) = (x − a)α + Φ(x)
Γ (α + 1)
onde a função Φ satisfaz
Φ(x) = O((x − a)µ+α )

quando x → a. Além disso,


 Hµ+α [a, b]
 se µ + α < 1
Φ∈ H ∗ [a, b] se µ + α = 1

H1 [a, b] se µ + α > 1

♡

Teorema 4.4
Seja α > 0, p > max{1, 1/α} e seja f ∈ Lp [a, b]. Então
Iaα f (x) = O((x − a)α−1/p )
quando x → a+ . Se α − 1/p ∈
/ N então Iaα f ∈ C [α−1/p] [a, b],e D[ α − 1/p]Iaα f ∈
Hα−1/p−⌊α−1/p⌋ [a, b]
♡

As demonstrações podem ser encontradas em [45].

O intercambio da operação de limite com o operador integral é um resultado conhecido, quando
α ∈ N. No próximo teorema discutimos a versão para o operador integral não inteiro.
Teorema 4.5
Sejam (fk )∞ k=1 uma sequência uniformemente convergente de funções contı́nuas em [a, b]. Então
é válida a seguinte igualdade.

Iaα lim fk (x) = lim Iaα fk (x)
k→∞ k→∞
Em particular, a sequência (Ia fk )∞
α
k=1 é uniformemente convergente.
♡

Demonstração Vamos considerar que f seja o limite da sequência (fk )∞ k=1 . O A primeira afirmação
segue diretamente do fato de que f é contı́nua. Para a segunda afirmação temos que se α = 0, o
operador Ia0 é o operador identidade e o resultado segue pelo fato da convergência uniforme da série.

40
4.1 Operador de Riemann-Liouville

Para α > 0 segue que

Z x
1
|Iaα fk (x) − Iaα f (x)| = |fk (t) − f (t)|(x − t)α−1 dt
Γ (α) a
Z b
1
≤ ∥fk − f ∥∞ (x − t)α−1 dt
Γ (α) a
Z b−a
1
≤ ∥fk − f ∥∞ (u)α−1 du
Γ (α) 0
1
≤ ∥fk − f ∥∞ (b − a)α
Γ (α + 1)

que converge uniformemente para zero, quando k → ∞, para todo x ∈ [a, b]

Com este resultado, podemos mostrar a conexão entre o operador integral não inteiro com a
derivação inteira de uma função analı́tica.
Corolário 4.1
Seja f função analı́tica em (a − h, a + h) para h > 0 e seja α > 0. Então
∞
X (−1)k (x − a)k+α k
Iaα f (x) = D f (x)
k=0
k!(α + k)Γ (α)
para a ≤ x < a + h/2, e
∞
X (x − a)k+α
Iaα f (x) = Dk f (a)
k=0
Γ (k + 1 + α)
para a ≤ x < a + h. Em particular, Iaα f é analı́tica em (a, a + h).
♡

Demonstração Sendo f uma função analı́tica, esta pode ser representada numa série de potências
sobre x. Além disso, x ∈ [a, a + h/2) e a série de potências converge em todo o intervalo de
integração. Segundo do Teorema 4.5 podemos trocar a ordem entre a soma e integração. Então,
usando a representação explicita para a integral não inteira a primeira parte é válida. A segunda parte
pode ser deduzida de forma análoga, expandindo a séria em a em vez de x. A analiticidade de Iaα f
segue imediatamente da segunda parte
Inspirado pelo resultado do Lema 2.2 podemos definir o operador diferencial não inteiro.
Definição 4.2 (Derivada Fracionária de Riemann-Liouville)
Seja α ∈ R+ e n = ⌈α⌉. O operador Daα definido por
n Z x
α n n−α 1 d
Da f (x) = D Ia f (x) = (x − t)n−α−1 f (t)dt
Γ (n − α) dx a
para a ≤ x ≤ b, é chamado de Riemann-Liouville de ordem α. Quando α = 0 temos o operador
identidade, ou seja, Da0 := I.
♣

Observe que, quando α ∈ N, o operador Daα coincide com o operador diferencial clássico. O
próximo lema é a versão não inteira do resultado do Lema 2.2.

41
4.1 Operador de Riemann-Liouville

Lema 4.1
Seja α ∈ R+ e seja n ∈ N tal n > α. Então,
Daα = Dn Ian−a
♡

Demonstração A suposição sobre n implica que n ≥ ⌈α⌉. Desta forma,

Dn Ian−α = D⌈α⌉ Dn−⌈α⌉ Ian−⌈α⌉ Ia⌈α⌉−α = D⌈α⌉ I ⌈α⌉−α = D⌈α⌉ Ia⌈α⌉−α = Daα
em razão das propriedades de semigrupo da integração não inteira e do fato que a diferenciação é a
inversa a esquerda do integração clássica
O operador diferencial não inteiro existe se certas condições são satisfeitas. No próximo lema
apresentamos essas condições.
Lema 4.2
Seja f ∈ A1 [a, b] e 0 < α < 1. Então Daα f existe em quase todo [a, b]. Além disso,Daα f ∈
Lp [a, b] para 1 ≤ p < 1/α e
Z x
α 1 f (a) ′ −α
Da f (x) = + f (t)(x − t) dt
Γ (1 − α) (x − a)α a
♡

Demonstração Ver Lema 2.12 (pag. 27) em [6]

A seguir são apresentados dois exemplos de aplicação do operador diferencial Riemann-Liouville.

Exemplo 4.1.1 Vamos considerar a função f (x) = (x − a)c , para algum c > −1 e α > 0. Desejamos
obter Daα f . Pela definição do operador temos.
Daα f (x) = D⌈α⌉ Ia⌈α⌉−α f (x)
Z x
⌈α⌉ 1 c ⌈α⌉−α−1
= D (t − a) (x − t) dt
Γ (α) a
Z 1
⌈α⌉ 1 ⌈α⌉−α+c β ⌈α⌉−α−1
= D (x − a) s (1 − s) ds
Γ (α) 0

⌈α⌉ Γ (c + 1) ⌈α⌉−α+c
= D (x − a)
Γ (⌈α⌉ − α + c + 1)
Γ (c + 1)
= D⌈α⌉ (x − a)⌈α⌉−α+c
Γ (⌈α⌉ − α + c + 1)

Segue-se que para o caso, onde (−α + c) ∈ N, o lado direito é simplesmente a ⌈α⌉-ésima derivada
de um polinômio de grau (⌈α⌉ − α + c) ∈ {0, 1, . . . , ⌈α⌉ − 1}, e desta forma a expressão se anula,
isto é,
Daα [(· − a)α−n ](x) = 0, para todo α > 0, n ∈ {1, 2, . . . , ⌈α⌉}

Quando (−α + c) ∈
/ N, segue que
Γ (c + 1)
Daα [(· − a)α−n ](x) = (x − a)c−α
Γ (c + 1 − α)

Observe que as duas últimas relações são simples generalizações do caso clássico. Mas também
revela um comportamento singular da derivada fracionária. No caso das funções constantes, o operador

42
4.1 Operador de Riemann-Liouville

Figura 4.1: Derivada de f (x) = 1, para α = 0.73.

diferencial de Riemann-Liouville apresenta solução não nula. Na Figura 4.1 temos a derivada de
f (x) = 1, para α = 0.73.
No próximo exemplo, veremos como a derivada não inteira, em geral, não pode ser obtida
diretamente do caso inteiro.

Exemplo 4.1.2 Seja f (x) = eλx , para algum λ > 0, e seja α > 0, α ∈
/ N. Então
∞
eλa X Γ (1 + k)Γ (1 − α) (λ(x − a))k
Daα f (x) = (x − a)−α .
Γ (1 − α) k=0
Γ (1 − α + k) k!
Para α = n ∈ N a expressão acima coincide com o resultado clássico.
Dan = eλx = λn eλx

Existe várias formas de definir o operador diferencial fracionário. Cada um com caracterı́sticas que
podem atender propósitos diferentes.
Devido a sua definição, o operador diferencial clássico satisfaz propriedade de semigrupo. No
teorema que segue veremos que o operador diferencial de Riemann-Liouville satisfaz propriedade
similar.
Teorema 4.6
Sejam α1 , α2 ≥ 0. Além disso, seja g ∈ L1 [a, b] e f = I α1 +α2 g. Então
Daα1 Daα2 f = Daα1 +α2 f.
♡

Demonstração Pela definição do operador diferencial de Riemann-Liouville, segue que,

Daα1 Daα2 f = Daα1 Daα2 Iaα1 +α2 g
= Da⌈α1 ⌉ Ia⌈α1 ⌉−α1 Da⌈α2 ⌉ Ia⌈α2 ⌉−α2 Iaα1 +α2 g
A propriedade de semigrupo do operador integral permite reescrever a expressão como
Daα1 Daα2 f = Da⌈α1 ⌉ Ia⌈α1 ⌉−α1 Da⌈α2 ⌉ I ⌈α2 ⌉+α1 g
= Da⌈α1 ⌉ Ia⌈α1 ⌉−α1 Da⌈α2 ⌉ Ia⌈α2 ⌉ Iaα1 g

43
4.1 Operador de Riemann-Liouville

O fato de que o operador diferencial clássico é o operador inverso a esquerda do operador integral e o
fato de que as ordens da integração e diferenciação envolvidas são números naturais segue que.
Daα1 Daα2 f = Da⌈α1 ⌉ Ia⌈α1 ⌉−α1 Iaα1 g = Da⌈α1 ⌉ Ia⌈α1 ⌉ g = g
A demonstração que Dα1 +α2 f = g segue de forma análoga
Aqui devemos fazer uma observação. O Teorema 4.4 impõe condições em f , sem as quais
não podemos garantir que a propriedade de semigrupo da derivada fracionária seja satisfeita. Nos
próximos exemplos,veremos casos em que, propriedades de semigrupos
Daα1 Daα2 f = Dα1 +α2 f (4.3)
Daα1 Daα2 f = D D fα2 α1
(4.4)
não são satisfeitas, quando f não satisfaz as condições do Teorema 4.4.
Exemplo 4.1 Seja f (x) = x−1/2 e α1 = α2 = 1/2. Pelo Exemplo 4.1.1 temos que se (⌈α⌉ − α + c) ∈
{0, 1, . . . , ⌈α⌉ − 1}, Daα [(· − a)α−n ](x) = 0. Assim, temos que D0α1 = D0α2 = 0. com isto
D0α1 D0α2 f (x) = 0, mas D0α1 +α2 f (x) = D01 f (x) = −(2x3/2 )−1 e desta forma (4.3) não é satisfeita.
√
Exemplo 4.2 Sendo f (x) = x1/2 e α1 = 1/2 e α2 = 3/2. Pelo Exemplo 4.1.1 temos que D0α1 = π/2
e D0α2 = 0. Desta forma D0α1 D0α2 f (x) = 0, mas D0α2 D0α1 f (x) = x−3/2 /4 e desta forma (4.4) não é
satisfeita.
Isto ocorre pelo fato de que f (x) pode estar no núcleo do operador diferencial Daα para algum
α ∈ R+ .
Versões análogas do Teorema 4.5 e do Corolário 4.1 para o operador diferencial são apresentadas
a seguir.
Teorema 4.7
Seja α > 0 e sejam (fk )∞
k=1 uma sequência uniformemente convergente de funções contı́nuas
em [a, b] e que Da fk existam para todo k. Além disso assuma que (Daα fk )∞
α
k=1 converge
uniformemente em [a + ϵ, b] para todo ϵ > 0. Então, para todo x ∈ (a, b], temos

α α
Da lim fk (x) = lim Da fk (x)
k→∞ k→∞
Em particular, a sequência (Iaα fk )∞
k=1 é uniformemente convergente. ♡

Com este resultado, podemos mostrar a conexão entre o operador diferencial não inteiro com a
derivação inteira de uma função analı́tica.
Corolário 4.2
Seja f função analı́tica em (a − h, a + h) para h > 0 e seja α > 0, α ∈/ N. Então
∞
X α (x − a)k−α

Daα f (x) = Dk f (x)
k=0
k Γ (k + 1 − α)
para a ≤ x < a + h/2, e
∞
X (x − a)k−α
Daα f (x) = Dk f (a)
k=0
Γ (k + 1 − α)
para a ≤ x < a + h. Em particular, Daα f é analı́tica em (a, a + h).
♡

44
4.1 Operador de Riemann-Liouville

O fato do operador diferencial não inteiro ser definido por Daα f = Dn I n−α f com n = ⌈α⌉, faz
com que a demonstração dos resultados acima seja feita de forma análoga feita nas demonstrações do
Teorema 4.5 e do Corolário 4.1.
Teorema 4.8
Sejam f1 e f2 funções definidas em [a, b] tal que Daα f1 e Daα f2 existam em quase todo [a, b].
Além disso, sejam c1 , c2 ∈ R. Então Daα (c1 f1 + c2 f2 ) existe em quase todo [a, b], e
Daα (c1 f1 + c2 f2 ) = c1 Daα f1 + c2 Daα f2
♡

Demonstração A demonstração decorre de forma imediada da definição de Daα que é linear

Dentre as regras de derivação, a regra generalizada da regra do produto para o operador diferencial
de Riemann-Liouville é apresentada no próximo teorema.
Teorema 4.9
Seja α > 0 e assuma que f e g sejam funções analı́ticas em (a − h, a + h). Então,
⌈α⌉ +∞
α
X α k α−k
X α
Da [f · g](x) = (D f )(x)(D g)(x) + (Dk f )(x)(Iak−α g)(x)
k=0
k k
k=⌈α⌉+1

para a < x < a + h/2

♡

Demonstração Do Corolário 4.2 temos que

+∞
α
X α (x − a)k−α
Da [f · g](x) = Dk [f · g](x)
k=0
k Γ (k + 1 − α)
para a ≤ x < a + h/2. Aplicamos a regra do produto clássica e trocamos a ordem dos somatórios,
isto segue que:
+∞ k
α
X α (x − a)k−α X k
Da [f · g](x) = Dj f (x)Dk−j g(x)
k=0
k Γ (k + 1 − α) j=0 j
+∞ X+∞
α (x − a)k−α k
X
= Dj f (x)Dk−j g(x)
j=0 k=j
k Γ (k + 1 − α) j
+∞ ∞
(x − a)ℓ+j−α

X
j
X α ℓ+j
= D f (x) Dℓ g(x)
j=0 ℓ=0
ℓ + j Γ (ℓ + j + 1 − α) j

Observe que, pela Definição 3.8

α ℓ+j α α−j
=
ℓ+j j j ℓ

45
4.1 Operador de Riemann-Liouville

Com isto segue que

+∞ ∞
(x − a)ℓ+j−α

X α X α−j
Daα [f · g](x) = j
D f (x) Dℓ g(x)
j=0
j ℓ=0
ℓ Γ (ℓ + j + 1 − α)
⌈α⌉ ∞
(x − a)ℓ+j−α

X α j
X α−j
= D f (x) Dℓ g(x)
j=0
j ℓ=0
ℓ Γ (ℓ + j + 1 − α)
+∞ ∞
(x − a)ℓ+j−α

X α j
X α−j
+ D f (x) Dℓ g(x)
j ℓ=0
ℓ Γ (ℓ + j + 1 − α)
j=⌈α⌉+1

Com os resultados dos Corolários 4.1 e 4.2 pode-se substituir os somatórios internos e chegar no
resultado desejado
Algumas observações sobre a diferenças entre a formulação clássica e a formulação fracionária
são necessárias. Observe que o Teorema 4.8 assume que ambas as funções são analı́ticas. Para a
função f isto se faz necessário para dar sentido a variação do k nos inteiros não negativos. Por outro
lado, para a função g apenas as derivadas até a ordem α são necessárias. Mas na demonstração,
assumimos que f · g é analı́tica e isto só é possı́vel se g for analı́tica. Como consequência duas
propriedades importantes da regra clássica não se transfere para o caso fracionário. Primeiro que a
fórmula fracionária de Leibniz não é simétrica devido ao segundo termo. Segundo, que necessitamos
de conhecer infinitamente mais derivadas da função f para o cálculo da derivada do produto f · g.
Enquanto no caso clássico é necessário conhecer apenas as derivadas até a n-ésima ordem. Outra
regra de derivação importante é a regra da cadeia. Neste caso a formulação fracionária tem uma
estrutura complexa [6], sem uso prático.
Outro ponto importante é a propriedade inversa entre os operadores de diferenciação e integração.
Teorema 4.10
Seja α > 0. Então, para toda f ∈ L1 [a, b], temos
Daα Iaα f = f
em quase todo lugar. Se existe uma função g ∈ L1 [a, b], tal que f = Iaα g, então
Iaα Daα f = f
em quase todo lugar. ♡

Demonstração Seja n = ⌈α⌉. Então pela definição de Daα e as propriedades de semigrupo da

integral fracionária e a inversa a esquerda do operador clássico,
Daα Iaα f (x) = Dn Ian−α Iaα f (x) = Dn Ian f (x) = f (x)
A segunda afirmação é uma consequência imediata do resultado acima. Sendo f = Iaα g, segue que
Iaα Daα f (x) = Iaα [Daα Iaα g] = Iaα g(x) = f (x)

A primeira afirmação mostra que o operador diferencial de Riemann-Liouville é de fato inverso à

esquerda do operador integral de Riemann-Liouville. Por outro lado o inverso à direita é obtido com

46
4.2 Operador de Caputo

a condição f = Iaα g. Se f não satisfizer esta condição, uma representação diferente para Iaα Daα f pode
ser obtida.
Teorema 4.11
Seja α > 0 e n = ⌈α⌉. Assumimos que f tal que Iaα f ∈ [a, b]. Então,
n−1
α α
X (x − a)α−k−1
Ia Da f (x) = f (x) − lim Dn−k−1 Ian−α f (z).
k=0
Γ (α − k) z→a+

Especificamente, para 0 < α < 1 segue que

(x − a)α−1
Iaα Daα f (x) − lim I 1−α f (z).
Γ (α) z→a+ a
♡

Com o resultado deste teorema segue a demostração da versão da expansão de Taylor no caso
fracionário.
Teorema 4.12
Seja α > 0 e n = ⌈α⌉. Assumimos que Ian−α f ∈ An [a, b]. Então
n−1
(x − a)α−n n−α
X (x − a)k+α−n
f (x) = lim+ Ia f (z)+ lim Dk+α−n f (z)+Iaα Daα f (x)
Γ (α − n + 1) z→a k=1
Γ (k + α − n + 1) z→a+
♡

Demonstração Através do Teorema 4.11 e fazendo a mudança de variável

n−1
X (x − a)α−k−1
f (x) = + Iaα Daα f (x)
k=0
Γ (α − k)
Fazendo k = 0, no somatório, obtemos o primeiro termo e ao restante aplicamos o resultado do Lema
4.1 para obter o resultado do teorema
O próximo resultado mostra a conexão da diferencial fracionária com a Transformada de Laplace.
Teorema 4.13
Sendo α > 0 e n = ⌈α⌉, seque que
n−1
X
(LD0α f )(s) α
= s (Lf )(s) − sn−k−1 Dk (I0n−α f )(0) (4.5)
k=0
♡

Observação O operador diferencial da Definição (4.1) também é chamado de operador diferencial de
Riemann-Liouville a esquerda, onde o operador integral tem como intervalo de integração [a, x], Neste
α
caso denotamos o operador diferencial por D[a,x] . Em contrapartida temos o operador diferencial a
α
direita D[x,a] , que na sua definição o operador integral tem como intervalo de integração [x, a], isto é
n Z a
α n a 1 d
D[x,a] f (x) = D In−α f (x) = (x − t)n−α−1 f (t)dt (4.6)
Γ (n − α) dx x
A escolha entre a derivação à direita ou à esquerda depende do contexto do problema e do que se
pretende calcular.

47
4.2 Operador de Caputo

4.2 Operador de Caputo

Uma nova definição para o operador diferencial não inteiro foi apresentado no trabalho de M.
Caputo [2], que ficou conhecido como operador diferencial de Caputo. Nesta seção iremos apresentar
algumas propriedades desse operador, iniciando com a sua definição.
Definição 4.3
Seja α ∈ R+ e seja n = ⌈α⌉. O operador C Daα definido por
Z x n
C α n−α n 1 n−α−1 d
Da f (x) = Ia D f (x) = (x − t) f (t)dt (4.7)
Γ (n − α) a dt
para a ≤ x ≤ b é chamado de operador diferencial de Caputo de ordem α. ♣

Observe que na definição do operador diferencial de Caputo, temos a inversão da da ordem que a
diferenciação e integração são aplicadas no operador diferencial de Riemann-Liouville. Este fato tem
um impacto bastante importante na estrutura da diferenciação fracionária, como pode ser verificado
com o exemplo que se segue:

Exemplo 4.2.1 Seja α ≥ 0, n = ⌈α⌉ e f (x) = (x − a)c para algum c ≥ 0. Então

(
C 0 se c ∈ {0, 1, . . . , n − 1}
Daα f (x) = Γ (c+1) c−α
Γ (c+1−α)
(x − a) se c > n − 1

Ao comparar com o Exemplo 4.1.1, notamos que o os operadores tem diferentes kernel, e que o
domı́nio dos dois operadores são diferentes. No entanto, um primeiro resultado conectando os dois
operadores é apresentado no próximo teorema.
Teorema 4.14
Seja α ≥ 0 e n = ⌈α⌉. Além disso, assuma que Daα f e f possui (n − 1) derivadas em a. Então
C
Daα f = Daα [f − Tn−1 [f ; a]]
em quase todo lugar, onde Tn−1 [f ; a] denota o Polinômio de Taylor de grau n − 1 para a função
f , centrado em a. ♡

Demonstração Do Lema 4.1 e da Definição 4.1 segue:

Daα [f − Tn−1 [f ; a]] = Dn Ian−α [f − Tn−1 [f ; a]]
Z x
dn (x − t)n−α−1
= (f (t) − Tn−1 [f ; a](t)) dt. (4.8)
dxn a Γ (n − α)
Aplicando a integração por partes, com dv = (x − t)n−α−1 dt e u = f (t) − Tn−1 [f ; a](t), segue que
Rx 1
a Γ (n−α)
(f (t) − Tn−1 [f ; a](t)) (x − t)n−α−1 dt.
1 t=a
= Γ n−α+1
[(f (t) − Tn−1 [f ; a](t))(x − t)n−α ]t=a
1
Rx
+ Γ (n−α+1) a
(Df (t) − DTn−1 [f ; a](t)(x − t)n−α dt
O termo que não depende da integral se anula, o primeiro fator quanto t = a e o segundo fator quando

48
4.2 Operador de Caputo

t = x. Desta forma
Ian−α [f − Tn−1 [f ; a]] = Ian−α+1 D[f − Tn−1 [f ; a]]

Com as suposições feitas e repetindo o processo n vezes, obtemos

Ian−α [f − Tn−1 [f ; a]] = Ia2n−α Dn [f − Tn−1 [f ; a]] (4.9)
= Ian Ian−α+1 Dn [f − Tn−1 [f ; a]]. (4.10)
Sendo Tn−1 [f ; a] um polinômio de grau n − 1, tem-se que Dn Tn−1 [f ; a] ≡ 0, assim sendo,
Ian−α [f − Tn−1 [f ; a]] = Ian Ian−α+1 Dn f.
Combinando com 4.8 obtém-se
Daα [f − Tn−1 [f ; a]] = Dn Ian Ian−α Dn f = Ian−α Dn f = C Daα f

Este teorema estabelece que o operador diferencial de Caputo de uma função f está definida se
o operador diferencial de Riemann-Liouville existe e se f é (n − 1)-vezes diferenciável no sentido
clássico, que garante a existência do polinômio de Taylor. Outra forma de expressar a relação entre os
operadores é dado pelo lema que se segue.
Lema 4.3
Seja α ≥ 0 e n = ⌈α⌉. Assuma que f é, tal que, ambos C Daα f e Daα existam. Então,
n−1
C α α
X Dk f (a)
Da f (x) = Da f (x) − (x − a)k−a
k=0
Γ (k − α + 1)
♡

Demonstração Segundo a definição do operador diferencial de Caputo e do Exemplo 4.1.1,

n−1
C α
X Dk f (a) α
Da f (x) = Daα f (x) − Da [(· − a)k ](x)
k=0
k!
n−1
X Dk f (a)
= Daα f (x) − (x − a)k−α
k=0
Γ (k − α + 1)

Uma consequência imediata desse resultado é apresentado no lema que segue.

Lema 4.4
Seja α ≥ 0 e n⌈α⌉. Assume que f é tal que que os operadores C Daα f e Daα f existam. Além
disso, seja Dk f (a) = 0, k = 0, 1, . . . , n − 1. Então,
C
Daα f = Daα f
♡

Isto é especialmente importante na solução de equações diferenciais de ordem não inteira. Se as

condições iniciais são homogêneas, as equações formuladas com o operador diferencial de Riemann-
Liouville coincide com a formulação com o operador diferencial de Caputo. Esta interação entre os
operadores determina que o operador Caputo é também inverso a esquerda do operador integral de
Riemann-Liouville.

49
4.2 Operador de Caputo

Teorema 4.15
Se f é contı́nua e α ≥ 0, então
C
Daα Iaα f = f
♡

Demonstração Sendo g = Iaα f . Pelo Teorema 4.3, segue que Dk g(a) = 0, para k = 0, 1, . . . , n − 1,
e desta forma
C α α
Da Ia f = C Daα g = Daα g = Daα g = Daα Iaα f = f

Em contrapartida o operador diferencial de Caputo não é inverso, a direita do operador integral

Riemann-Liouville.
Teorema 4.16
Seja α ≥, n = ⌈α⌉, e f ∈ An [a, b]. Então,
n−1
X Dk f (a)
Iaα C Daα f = f (x) − (x − a)a
k=0
k!
♡

Demonstração Sendo g = Iaα f . Pelo Teorema 4.3, segue que Dk g(a) = 0, para k = 0, 1, . . . , n − 1,
e desta forma
C α α
Da Ia f = C Daα g = Daα g = Daα g = Daα Iaα f = f

A versão da expansão de Taylor para o operador diferencial de Caputo é apresentada a seguir.

Corolário 4.3
Seja α > 0 e n = ⌈α⌉. Seja f é tal que f ∈ An [a, b]. Então
n−1
X Dk f (a)
f (x) = (x − a)k + Iaα C Daα f (x)
k=0
k!
♡

Comparando com o resultado do Teorema 4.12, que apresenta expansão de Taylor para o operador
diferencial Riemann-Liouville, observa-se que a expansão do operador Caputo tem uma estrutura bem
simples.
O resultados a seguir, são as versões das regras de derivação para o operador diferencial de
Caputo.
Teorema 4.17
Seja f1 e f2 funções definidas em [a, b] tal que C Daα f1 e C Daα f2 existem em quase todo lugar.
Seja c1 , c2 ∈ R. Então,
C
Daα (c1 f1 + c2 f2 ) = c1 C Daα f1 + c2 C Daα f2
♡

Demonstração A definição do operador garante a linearidade.

50
4.2 Operador de Caputo

Teorema 4.18 (Fórmula de Leibniz)

Seja 0 < α < 1 e assuma que f e g são analı́ticas em (a − h, a + h). Então,
∞
C α (x − a)−α C α
X α
Da [f · g] = g(a)(f (x) − f (a)) + ( Da g(x)f (x) + (Iak−α g(x))C Daα f (x)
Γ (1 − α) k=1
k
♡

Demonstração Aplicando a definição do operador diferencial de Caputo segue que

C
Daα [f · g] = Daα [f · g − f (a)g(a)] = Daα [f · g] − f (a)g(a)Daα [1]
Aplicamos a fórmula de Leibniz para o operador diferencial de Riemann-Liouville e obtemos
∞
C α α
X α
Da [f · g] = f (Da g) + (Dak f )(Iak−α g) − f (a)g(a)Daα [1]
k=1
k
Somando e subtraindo o termo f (a)g(a)(Daα [1] e reagrupando os termos segue que
∞
C α α
X α
Da [f · g] = f (Da [g − g(a)]) + (Dak f )(Iak−α g) + g(a)(f − f (a))Daα [1]
k=1
k
∞
X α C α
= f (Daα g) + ( Da f )(Iak−α g) + g(a)(f − f (a))Daα [1]
k=1
k

onde usamos o fato que, para k ∈ N, Dak = Dk = C Daα . Uma expressão para Daα [1] é obtida pelo
Exemplo 4.1.1, que completa a demonstração
Nos próximos dois lemas veremos outras diferenças significativas entre os operadores diferenciais
de Riemann-Liouville e de Caputo.
Lema 4.5
Seja α > 0, α ∈ / N e n = ⌈α⌉. Além disso, assuma que f ∈ C n [a, b]. Então C Daα f é contı́nua
em [a, b] e C Daα f (a) = 0
♡

Demonstração Pela definição e Teorema 4.14, C Daα f = Ian−α Dn f . O resultado segue pelo Teorema
4.3 porque Dn f é contı́nua por hipótese
As condições do lema anterior pode ser relaxadas.
Lema 4.6
Seja α > 0, α ∈ / N e n = ⌈α⌉. Além disso, assuma que f ∈ An [a, b] e assumimos que
C ē
Da f ∈ C[a, b] para algum ē ∈ [α, n). Então C Daα f é contı́nua em [a, b] e C Daα f (a) = 0
♡

Demonstração Pela definição e Teorema 4.14, e a propriedade de semigrupo do opearador integral

não inteiro,
C α ē
Da f = Ian−α Dn f = Iaē−α Ian−ē Dn f = Iaē−α C Da f

O resultado segue pelo Teorema 4.3

A Transformada de Laplace do operador diferencial de Caputo é dado por:

51
4.3 Operador de Grünvald-Letnikov

Teorema 4.19
Sendo α > 0 e n = ⌈α⌉, seque que
n−1
X
C
(L Daα f )(s) α
= s (Lf )(s) − sα−k−1 (Dk f )(0) (4.11)
k=0
em particular, se 0 < α ≤ 1, teremos
(LC Daα f )(s) = sα (Lf )(s) − sα−1 f (0)
♡

4.3 Operador de Grünvald-Letnikov

O operador diferencial não inteiro de Grünwald-Letnikov leva o nome dos matemáticos Augustin
Louis Cauchy Grünwald e Aleksandr Vasilevich Letnikov, que contribuı́ram independentemente para
o desenvolvimento dessa abordagem.
A história começa com Augustin Louis Cauchy Grünwald, um matemático austrı́aco, que traba-
lhou no século XIX. Ele fez importantes contribuições para o campo da análise matemática, incluindo
o estudo de séries e equações diferenciais. Em seu trabalho, Grünwald investigou a ideia de generalizar
a noção de derivada para ordens não inteiras.
Mais tarde, no final do século XIX, Aleksandr Vasilevich Letnikov, um matemático russo, também
começou a explorar a mesma questão de derivadas não inteiras. Letnikov desenvolveu uma abordagem
usando diferenças finitas para calcular as derivadas em pontos discretos de uma função.
No entanto, o trabalho de Grünwald e Letnikov permaneceu relativamente desconhecido até o
século XX, quando o campo da análise fracionária começou a ganhar mais atenção. Os estudos de
Grünwald e Letnikov foram redescobertos e seu trabalho passou a ser considerado fundamental para
o desenvolvimento dessa área.
Definição 4.4
Seja α ∈ R+ . O operador GL Daα , definido por
m
(∆αh f )(x)

GL α 1 X k α
Da f (x) = lim = lim α (−1) f (x − kh) onde mh = x − a (4.12)
h→0 hα h→0 h
k=0
k
para a ≤ x ≤ b, é chamado de operador diferencial não inteiro de Grünwald-Letnikov de ordem
α. ♣

Nesta definição o termo (∆αh f )(x) é a formulação não inteira do operador diferença. A definição
é válida para uma função arbitrária f (x), mas a convergência da soma infinita não é garantida para toda
função. No entanto, aplicando o operador diferencial de Grünwald-Letnikov em equações diferenciais
de ordem não inteira, podemos definir a solução y(x) formalmente de zero a menos infinito. Como
resultado. a soma infinita na Definição 4.4 se torna uma soma finita. Também é possı́vel definir o

52
4.3 Operador de Grünvald-Letnikov

operador integral de Grünwald-Letnikov não inteiro, trocando α por −α, obtendo

m
k −α
X
GL α α
Ia f (x) = lim h (−1) f (x − kh)
h→0
k=0
k
m
hα X Γ (k + α)
= lim f (x − kh) onde mh = x − a
h→0 Γ (α) Γ (k + 1)
k=0
O resultado do próximo teorema fornece uma visão da conexão entre a derivada de Grünwald-
Letnikov e as derivadas não inteiras definidas anteriormente.
Teorema 4.20
Seja α ≥ 0, n = ⌈α⌉ e f ∈ C n [a, b], então
n−1 (k) Z x
GL α
X f (a)(x − a)k−α 1
Da f (x) = + (x − t)n−1−α f (n) (t)dt (4.13)
k=0
Γ (k + 1 − α) Γ (n − α) a
♡

Corolário 4.4
Seja α ≥ 0, n = ⌈α⌉ e f ∈ C n [a, b], então
GL
Daα f (x) = Tn−1 [f ; a](x) + C Daα f (x) = Daα f (x).
♡

Demonstração O resultado da primeira igualdade segue do Teorema 4.20 e a segunda igualdade do

Teorema 4.14
A condição de que f ∈ C n [a, b] restringe a classe de funções para qual a conexão entre os
operadores possa ser estabelecida. Além disso, a definição do operador Caputo já exige a existência
da n-ésima derivada da função f , uma restrição não mais forte que f ∈ C n [a, b].
Para o desenvolvimento de algoritmos numéricos, é interessante saber o que ocorre com a soma
finita.
Teorema 4.21
Seja f ∈ C n [0, X], α ≥ 0, n = ⌈α⌉, e X/h = m ∈ N. Então o operador diferencial de
Grünvald-Letnikov finito, centrado em zero.
m
GL α 1 X k α
F Da f (x) = α (−1) f (x − kh)
h k=0 k
gera uma aproximação para o operador diferencial Riemann-Liouville Dα se e somente se
f (0) = x0 , isto é
GL α
F Da f (x) = Dα f (x) + O(h) + O(x0 )
♡

Demonstração O resultado segue imediatamente do Corolário 4.4

Para o nosso propósito de descrever os métodos numéricos, as três versões dos operadores
diferenciais apresentados serão suficientes, mas devemos ressaltar que existem outras maneiras de
definir operadores diferenciais não inteiro.

53
Capı́tulo 5 Equação Diferencial Ordinária
Fracionária
Neste capı́tulo iremos discutir a existência e unicidade das equações diferenciais ordinárias que
envolvem operadores diferenciais fracionário do tipo Riemann-Liouville e do tipo Caputo. A nossa
abordagem se restringe aos problemas de valor inicial (P.V.I.) e sem perda de generalidade assumimos
que os operadores são desenvolvidos no ponto x = 0. Desta forma na anotação dos operadores iremos
suprir o subı́ndice a.
Como no caso dos operadores diferenciais e integrais inteiros, esquemas numéricos, baseados
no resultados teóricos, tem sido uma ferramenta atrativa no campo da engenharia para descrever
diferentes tipos de modelos. Iremos focar nos esquemas numéricos para equações diferenciais ordinária
fracionária de cada tipo apresentado no Capı́tulo 4.

5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

Iniciamos com a definição de uma Equação Diferencial Ordinária Fracionária (EDOF ) do tipo
Riemann-Liouville e do tipo Caputo.
Definição 5.1
/ N, n = ⌈α⌉ e f : A ⊆ R2 → R. Então
Seja α > 0, α ∈
Dα y(x) = f (x, y(x)). (5.1)
é chamada de equação diferencial fracionária ordinária do tipo Riemann-Liouville. Como
condições iniciais para este tipo de EDOF nós usamos
Dα−k y(x) = bk , (k = 1, 2, . . . , n − 1), lim+ I n−α y(z) = bn . (5.2)
z→0
De forma análoga temos
C
Dα y(x) = f (x, y(x)). (5.3)

é chamada de equação diferencial fracionária ordinária do tipo Caputo. Neste caso, as

condições iniciais para este tipo de EDOF são dadas por
Dk y(x) = bk , (k = 1, 2, . . . , n − 1) (5.4)
♣

O uso de condições iniciais diferentes é necessários para garantir a unicidade da solução, a qual
será demonstrada no próximo teorema.
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

Teorema 5.1
/ N e n = ⌈α⌉. Além disso, seja K > 0, h∗ > 0 e b1 , b2 , . . . , bm ∈ R. Defini-se
Seja α > 0, α ∈
G := {(x, y) ∈ R2 : 0 ≤ x ≤ h∗ , y ∈ R para x = 0 e
Xn
|xn−α y − bk xn−k /Γ (α − k + 1)| < K},
k=1
e assuma que a função f : G → R é contı́nua e limitada em G e que esta satisfaz a condição
de Lipschitz com respeito a segunda variável, isto é, existe uma constante L > 0 tal que, para
todo (x, y1 ), (x, y2 ) ∈ G, temos
|f (x, y1 ) − f (x, y2 )| < L|y1 − y2 |.
Então a EDOF do tipo Riemann-Liouville (5.1) com as condições iniciais (5.2) tem uma única
solução contı́nua y ∈ C(0, h], onde h := min{h∗ , h1 , h2 }, onde
1/n
Γ (2α − n + 1) Γ (α + 1)K
h1 < , h2 = .
(Γ (α − n + 1)L)1/α M
♡

Para a equação diferencial do tipo Caputo, temos um resultado similar.

Teorema 5.2
/ N e n = ⌈α⌉. Além disso, seja K > 0, h∗ > 0 e b0 , b2 , . . . , bn−1 ∈ R. Defini-se
Seja α > 0, α ∈
G := [0, h∗ ] × [b0 − K, b0 + K],
e assuma que a função f : G → R é contı́nua. Então, existe um h > 0 e a função y ∈ C[0, h]
é solução da equação diferencial não inteira do tipo Caputo (5.3) com as condições iniciais
(5.4). Para o caso, α ∈ (0, 1) o parâmetro h é dado por
h := min{h∗ , (Γ (α + 1)KM )1/n }, com M := sup |f (x, z)|.
(x,z)∈G

G e que esta satisfaz a condição de Lipschitz com respeito a segunda variável, isto é, existe uma
constante L > 0 tal que, para todo (x, y1 ), (x, y2 ) ∈ G, temos
♡

Os resultados apresentados se assemelham aos resultados equivalentes no caso clássico. A

demonstração é feita de forma semelhante. Não faremos a demonstração, mas iremos apresentar, que
ambos os tipos, podem ser formuladas como equações integrais, ou seja como equações integrais de
Volterra.
Teorema 5.3
1. Sobre as hipóteses do Teorema 5.1 a função y ∈ C(0, h] é solução da equação diferencial
de Riemann-Liouville (5.1) com as condições iniciais (5.2), se e somente se é solução da
equação integral de Volterra de segunda ordem.
n Z x
X bk xα−k 1
y(x) = + (x − t)α−1 f (t, y(t))dt. (5.5)
k=1
Γ (α − k + 1) Γ (α) 0

2. Sobre as hipóteses do Teorema 5.2 a função y ∈ C[0, h] é solução da equação diferencial

de Caputo (5.3) com as condições iniciais (5.4), se e somente se é solução da equação

55
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

integral de Volterra de segunda ordem.

n−1 Z x
X bk x k 1
y(x) = + (x − t)α−1 f (t, y(t))dt. (5.6)
k=1
k! Γ (α) 0
♡

Observação Antes de fazermos a demostração, algumas observações sobre os diferentes tipos de
equações devem ser feitas.
a) No Teorema 5.1 assumimos que y é contı́nua em (0, h] e não em [0, h] como no caso clássico.
A equação (5.5) justifica este fato. Se assumirmos que y é contı́nua em [0, h] então o lado
esquerdo da equação integral é continua neste intervalo, assim como o lado direito da equação
integral (a continuidade de f justifica esse fato). Desta forma a soma será contı́nua em [0, h]
também. Em vista da definição de n, podemos observar, que de fato, as somas são contı́nuas
em [0, h], para k = 1, 2, . . . , n − 1, mas para k = n é ilimitado quando x → 0, porque n > α,
se α ∈/ N, a menos que bn = 0. Desta forma, y não é contı́nua na origem, a menos que bn = 0.
b) Com argumentos semelhantes para a equação integral (5.6), podemos observar que y é contı́nua
no intervalo [0, h] no Teorema 5.2: A soma no lado direito da equação é contı́nua, uma vez
que temos monômios de ordem inteira, não negativa, ao contrário do que ocorre em (5.5) com
k = n, onde o termo xα−n . é único em zero.
c) Comparando as condições iniciais (5.2), para as equações de Riemann-Liouville, com as
condições iniciais (5.4), para as equações de Caputo, observa-se no caso de Caputo, estas
condições solução formuladas como no caso clássico. Embora, em ambos os casos, a existência
e unicidade são garantidas, no caso de Caputo temos um impacto importante nas aplicações.
Em geral as condições iniciais para o tipo de Riemann-Liouville, não tem significado óbvio e
não podem ser medidos.
A seguir iremos fazer a Demonstração do Teorema 5.3.
Demonstração (Teorema 5.3) 1. Para a primeira parte do teorema, vamos assumir que y é uma
solução da equação integral. Reescrevendo a equação, usando a notação do operador integral temos
n
X bk xα−k
y(x) = + I0α f (·, y(·))(x).
k=1
Γ (α − k + 1)
Aplicamos o operador diferencial D0α em ambos os lados, verificamos que y também é solução da
equação diferencial. (Consequência do Exemplo 4.1.1 e Teorema 4.10). Com relação as condições
iniciais, inicialmente consideremos os casos em que 1 ≤ k ≤ n − 1. Aplicando o operador D0α−k na
equação de Volterra, obtén-se
n
α−k
X bj D0α−k (·)α−j (x)
D0 y(x) = + D0α−k I0α−k I0k f (·, y(·))(x).
j=1
Γ (α − j + 1)
em vista da propriedade de semigrupo da integração não inteira. Pelo Exemplo 4.1.1 conclui-se que
os termos do somatório se anulam para j > k. Além disso, pelo mesmo exemplo, os termos do
somatório, para j < k, se anulam quando x = 0, desta forma
bk D0α−k (·)α−k (x)
D0α−k y(x) = + I0k f (·, y(·))(x).
Γ (α − k + 1)

56
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

Como a integral se anula para k ≥ 1 e pelo Exemplo 4.1.1, obtemos que D0α−k (·)α−k (x) = Γ (α−k+1),
resultando que D0α−k y(0) = bk , ou seja satisfaz as condições iniciais para 1 ≤ k ≤ n − 1.
Para k = n, aplicamos o operador I0n em ambos os lados da equação integral e aplicando o limite
com z → 0, todos os termos do somatório se anulam, com exceção do termo de ı́ndice n. A integral
I0n−α I0α f (·, y(·))(z) = I0n f (·, y(·))(z) também se anula quando z → 0. Com isto, segue que
bn I n−α (·)α−n (z)
lim+ I0n y(z) = lim z → 0+I0n−α 0 = bn .
z→0 Γ (α − n + 1)
Desta forma y é solução do P.V.I.
Se y é uma solução contı́nua do P.V.I. então definimos z(x) := f (x, y(x)). Supondo que z seja
contı́nua e z(x) = f (x, y(x)) = Dα y(x) = Dn I0n−α y(x). Desta forma Dn I0n−α y também é contı́nua.,
isto é, I0n−α y ∈ C n (0, h]. Isto permite aplicar o Teorema 4.11, obtendo
Xn X n
α α α−k α
y(x) = I0 D0 y(x) + ck x = I0 f (·, y(·))(x) + ck xα−k ,
k=1 k=1
onde ck são determinados pelas condições iniciais, sendo
bk
ck = .
Γ (α − k + 1)
2. Pela definição do operador diferencial de Caputo, podemos reescrever a equação diferencial
fracionária (5.1) como
f (x, y(x)) = C Daα = Dα (y − Tn−1 [y, 0])(x) = Dn I n−α (y − Tn−1 [y; 0])(x)
No nosso contexto, as funções são contı́nuas e podemos aplicar o operador de integração n-vezes em
ambos dos lados, obtendo
I n f (x, y(x)) = I 1−α = (y(x) − Tn−1 [y, 0])(x) + q(x),
onde q(x) é um polinômio de grau menor ou igual a n − 1. Como f (x, y(x)) é contı́nua, I n f (x, y(x))
tem um zero de multiplicidade n na origem. Pela definição do polinômio de Taylor isto também é
verdade para o termo I 1−α = (y(x) − Tn−1 [y, 0])(x) e também devido a equação para q(x). Assim
q(x) = 0 e consequentemente
I n f (x, y(x)) = I 1−α = (y(x) − Tn−1 [y, 0])(x).
Aplicando o operador diferencial de Riemann-Liouville Dn−α em ambos os lados da equação segue
que
I α f (x, y(x)) = y(x) − Tn−1 [y, 0](x),

que é a equação integral de Volterra (5.6). A outra direção da igualdade pode ser provada pela aplicação
do operador diferencial C Daα em ambos os lados da equação integral de Volterra (5.6) resultando em
n−1 k
! Z x
C α C α
X x C α 1
Da y(x) = Da bk + Da (x − t)α−1 f (t, y(t))dt.
k=0
k! Γ (α) 0

Como o polinômio está dentro do kernel do operador C Daα apenas o termo C Daα I α f (x, y(x)) permanece
do lado direito, de modo que a equação se reduz a
C
Daα y(x) = f (x, y(x)),
que é a equação diferencial fracionária (5.3). Para verificar as condições iniciais, derivamos a equação

57
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

(5.6) k-vezes, com k = 0, 1, . . . , n − 1 usando a regra de Leibniz para a diferenciação de uma integral,
dando
n−k−1
X xj−k c(k) x
Z
(k)
y (x) = bk + (x − t)α−k−1 f (t, y(t))dt,
j=k
(j − k)! Γ (α) 0

onde c(k) é uma constante formada por termos que surgem da derivação, onde
Yk
c(k) = (α − k).
j=1
Para os casos k < n − 1 apenas o termo bk não se anula para x = 0, porque os outros termos da soma
são nulos e a integração de 0 a 0 é nula. No caso k = n − 1 a integral contém singularidade na origem
e desta forma justificamos o resultado pela integral imprópria.
Z x
lim (x − t)α−k−1 f (t, y(t))dt = 0
x→0 0
e desta forma, apenas o termo bn−1 não se anula. Com isto conclui-se que as condições iniciais são
válidas
Este resultado permite provar resultados análogos aos apresentados pelos Teoremas 5.1 e 5.2,
que estabelece os resultados em termos da correspondente equação integral.
Lema 5.1
Sobre as condições do Teorema 5.1 a equação de Volterra
n Z x
X bk xα−k 1
y(x) = + (x − t)α−1 f (t, y(t))dt (5.7)
k=1
Γ (α − k + 1) Γ (α) 0

possui uma única solução y ∈ C(0, h].

♡

Demonstração Iniciamos com a definição do conjunto

( n
)
X bk xα−k
B := y ∈ C(0, h] : sup xn−α y(x) − ≤K
0≤x≤h
k=1
Γ (α − k + 1)
e neste conjunto definimos o seguinte operador
n Z x
X bk xα−k 1
Ay(x) = + (x − t)α−1 f (t, y(t))dt
k=1
Γ (α − k + 1) Γ (α) 0
Desta forma, para y ∈ B, Ay é uma função contı́nua em (0, h]. Além disso,
n
bk xα−k xn−α x
X Z
n−α
x Ay(x) − = (x − t)α−1 f (f, y(t))dt
k=1
Γ (α − k + 1) Γ (α) 0
n−α Z x
x
≤ M (x − t)α−1 dt
Γ (α) 0
xn−α xα
≤ M ≤K
Γ (α) a

para x ∈ (0, h], onde a última desigualdade se deve da definição de h, Isso mostra que Ay ∈ B, se
y ∈ B, isto é, o operador A de B em B.

58
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

Consideremos o conjunto

n−α
B̂ := y ∈ C(0, h] : sup x y(x) ≤ ∞
0≤x≤h

e neste conjunto definimos a norma ∥ · ∥B̂ por

∥y∥B̂ := sup |xn−α y(x)|.
0<x≤h

Prova-se que B̂, com essa norma, é um espaço linear normado e B é um subespaço de B̂.
A equação de Volterra, pode ser reescrita, de forma reduzida como
y = Ay.
Consequentemente, temos que mostrar que o operador A tem um único ponto fixo. Para este propósito,
usamos o Teorema 2.1. Neste contexto, provaremos que para y, ỹ ∈ B,
α j
j j Lh Γ (α − n + 1)
∥A y − A ỹ∥B̂ ≤ ∥y − ỹ∥B̂ (5.8)
Γ (2α − n + 1)
A demonstração segue pelo processo de indução finita. Para j = 0 segue que
∥A0 y − A0 ỹ∥B̂ = ∥y − ỹ∥B̂
Considerando válido para j − 1, temos
∥Aj y − Aj ỹ∥B̂ = sup xn−α (Aj y(x) − Aj ỹ(x)
0<x≤h

= sup xn−α (AAj−1 y(x) − AAj−1 ỹ(x)

0<x≤h
Z x
xn−α
(x − t)α−1 f (t, Aj−1 y(t)) − f (t, Aj−1 ỹ(t)) dt

= sup
0<x≤h Γ (α) 0
xn−α x
Z
≤ sup (x − t)α−1 f (t, Aj−1 y(t)) − f (t, Aj−1 ỹ(t)) dt
0<x≤h Γ (α) 0
n−α Z x
x
≤ sup (x − t)α−1 Aj−1 y(t) − Aj−1 ỹ(t) dt
0<x≤h Γ (α) 0

pela definição do operador A e a condição de Lipschitz da função f . Continuando a estimativa

Z x
j j L n−α
∥A y − A ỹ∥B̂ ≤ sup x (x − t)α−1 Aj−1 y(t) − Aj−1 ỹ(t) dt
Γ (α) 0<x≤h
Z0 x
L
≤ sup xn−α (x − t)α−1 tα−n tn−α Aj−1 y(t) − Aj−1 ỹ(t) dt
Γ (α) 0<x≤h 0
Z x
L j−1 j−1 n−α
≤ A y(t) − A ỹ(t) B̂ sup x (x − t)α−1 tα−n dt
Γ (α) 0<x≤h 0
L Γ (α)Γ (α − n + 1) α
= Aj−1 y(t) − Aj−1 ỹ(t) B̂ sup x
Γ (α) 0<x≤h Γ (2α − n + 1)
Lhα Γ (α − n + 1)
= Aj−1 y(t) − Aj−1 ỹ(t) B̂
Γ (2α − n + 1)

E usando a hipótese de indução provamos para j. Aplicando o resultado do Teorema 2.1, com αj = γj ,

59
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

onde
Lhα Γ (α − n + 1)
γ= .
Γ (2α − n + 1)
∞
X
Falta demonstrar que que a série αj converge. Mas isto decorre de que h ≤ h̃ e pela definição de
j=0
h̃ segue que γ < 1. Pelo teorema do ponto fixo seque que a nossa equação integral tem uma solução
única
Lema 5.2
Sobre as condições do Teorema 5.2 a equação de Volterra
n−1 k Z x
X x k 1
y(x) = D y(0) + (x − t)α−1 f (t, y(t))dt (5.9)
k=0
k! Γ (α) 0

possui uma única solução y ∈ C[0, h].

♡

Demonstração A demonstração é dividida em duas partes. Na primeira, consideramos o caso de

α > 1 e na segunda o caso em que α ∈ (0, 1). A razão disto, é que a equação de Volterra possui uma
singularidade no kernel (x − t)α−1 no caso α ∈ (0, 1) e é contı́nuo nos outros casos.
Consideremos α > 1. A equação (5.9) possuı́ um kernel contı́nuo resultando numa função
contı́nua na solução da integral. Assim a existência da solução segue os métodos padrão das equações
de Volterra [31].
No caso α ∈ (0, 1), a equação (5.9) se reduz a
Z x
1
y(x) = b0 + (x − t)α−1 f (f, y(t))dt. (5.10)
Γ (α) 0
Para a demonstração de existência de solução nós introduzimos o conjunto
U := {y ∈ C[0, h] : ∥y − b0 ∥∞ ≤ K}
O conjunto U é um subespaço de Banach convexo e fechado de todas as funções contı́nuas em [0, h],
equipado com a norma de Chebyshev. Desde que a função constante y = b0 pertence a U . Em U
definimos o operador
Z x
1
(Ay)(x) := b0 + (x − t)α−1 f (t, y(t))dt. (5.11)
Γ (α) 0
Com isto reescrevemos a equação (5.10) com sendo um problema de ponto fixo
y = Ay
Desta forma, vamos investigar as propriedades do operador A.

60
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

Note que, para 0 ≤ x1 ≤ x2 ≤ h,

Z x1 Z x2
1 α−1
∥(Ay)(x1 ) − (Ay)(x2 )∥ = (x1 − t) f (t, y(t))dt − (x2 − t)α−1 f (t, y(t))dt
Γ (α) 0 0
Z x1 Z x2
1 α−1 α−1
= ((x1 − t) − (x2 − t) )f (t, y(t))dt + (x2 − t)α−1 f (t, y(t))dt
Γ (α) 0 x1
Z x1 Z x2
∥f ∥∞ α−1 α−1
≤ ((x1 − t) − (x2 − t) ) + (x2 − t)α−1 dt
Γ (α) 0 x1
∥f ∥∞
= (2(x2 − x1 )α + xα1 − xα2 ),
Γ (α + 1)
provando que Ay é uma função contı́nua. Além disso, para y ∈ U e x ∈ [0, h], temos
Z x
(0) 1
|(Ay)(x) − y0 | = (x − t)α−1 f (t, y(t))dt
Γ (α) 0
1
≤ ∥f ∥∞ xα
Γ (α + 1)
1
≤ ∥f ∥∞ hα
Γ (α + 1)
1 KΓ (α + 1)
≤ ∥f ∥∞ =K
Γ (α + 1) ∥f ∥∞

Desta forma, mostramos que para y ∈ U , Ay ∈ U , isto é, operador A mapeia U em si mesmo.
Desejamos aplicar o Teorema 2.6 (Teorema do ponto fixo de Schauder) e para isto, necessitamos
mostrar que A(U ) é relativamente compacto. Isto pode ser feito por intermédio do Teorema 2.7
(Teorema de Arzelà-Ascoli). Para z ∈ A(U ), segue que, para todo x ∈ [0, h],
|z| = |(Ay)(x)|
Z x
1
≤ |b0 | + (x − t)α−1 |f (t, y(t))|dt
Γ (α) 0
1
≤ |b0 | + ∥f ∥∞ hα ,
Γ (α + 1)
o que demonstra ser limitado. Além disso, para 0 ≤ x1 ≤ x2 ≤ h e se |x2 − x1 | < δ, então
∥f ∥∞ α
|(Ay)(x1 ) − (Ay)(x2 )| ≤ 2 δ
Γ (α + 1)
Observe que a expressão a direita é independente de y, mostrando que A(U ) é equicontinuo. Então o
Teorema de Arzelà-Ascoli garante que A(U ) é relativamente compacto. Desta forma o Teorema do
ponto fixo de Schauder garante que A tem um ponto fixo. Pela definição de A, este ponto fixo é a
solução do p.v.i. Basta provar a unicidade do ponto fixo.
Para isto, vamos provar que, para todo j ∈ N0 e para todo x ∈ [0, h] temos
(Lxα )j
∥Aj y − Aj ỹ∥L∞ [0,x] ≤ ∥y − ỹ∥L∞ [0,x] (5.12)
Γ (1 + αj)
Isto pode ser visto por indução. Para j = 0 o resultado é trivial. Supondo verdadeiro para j − 1 segue.

61
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

∥Aj y − Aj ỹ∥L∞ [0,x] = ∥A(Aj−1 y) − A(Aj−1 ỹ)∥L∞ [0,x]

Z w
1
(w − t)α−1 f (t, Aj−1 y(t)) − f (t, Aj−1 ỹ(t)) dt

= sup
Γ (α) 0≤w≤x 0

No próximo passo, usamos a suposição de Lipschitz em f e a hipótese de indução, obtendo

Z w
j j L
∥A y − A ỹ∥L∞ [0,x] ≤ sup (w − t)α−1 Aj−1 y(t) − Aj−1 ỹ(t) dt
Γ (α) 0≤w≤x 0
Z x
L
≤ (x − t)α−1 sup Aj−1 y(w) − Aj−1 ỹ(w) dt
Γ (α) 0 0≤w≤t
j Z x
L
≤ (x − t)α−1 tα(j−1) sup |y(w) − ỹ(w)| dt
Γ (α)Γ (1 + α(j − 1)) 0 0≤w≤t
j Z x
L
≤ sup |y(w) − ỹ(w)| (x − t)α−1 tα(j−1) dt
Γ (α)Γ (1 + α(j − 1)) 0≤w≤x 0
j
L Γ (α)Γ (1 + α(j − 1)) αj
= ∥y(w) − ỹ(w)∥L∞ [0,x] x
Γ (α)Γ (1 + α(j − 1)) Γ (1 + αj)

o que resulta na equação (5.12). Aplicando a norma de Chebyshev no intervalo fundamental [0, h],

(Lxα )j
∥Aj y − Aj ỹ∥∞ ≤ ∥y − ỹ∥∞
Γ (1 + αj)
Com isto mostramos que o operador A satisfaz as condições do Teorema 2.1, com
(Lxα )j
αj =
Γ (1 + αj)
Para aplicar o teorema, é necessário verificar que a série
X∞
S= αj
j=0
converge. Isto, no entanto, segue imediatamente da propriedade 3, do Teorema 3.4. Com isto podemos
aplicar o Teorema do ponto fixo de Weissinger e deduzir a unicidade da solução da equação diferencial

A série S da demonstração do Lema 5.2 nada mais é do que a função de Mittag-Leffler, Eα (Lhα ).
A função de Mittag-Leffler não apenas surge neste contexto, ela realmente tem um papel fundamental
no campo do Cálculo Fracionário, como exemplificado no próximo teorema.
Teorema 5.4
Seja α > 0, n = ⌈α⌉ e λ ∈ R. A solução do P.V.I.
C
Daα y(x) = λy(x), y(0) = 1y k (0) = 0(k = 1, 2, . . . , n − 1)
é dada por
y(x) = Eα (λxα ), x ≥ 0

Em outras palavras, as autofunções do operador diferencial de Caputo pode ser expressas em

62
5.1 Equações do tipo Riemann-Liouville e do tipo Caputo

termos das funções de Mittag-Leffler.

♡

Demonstração Em resultado anterior, foi demonstrado a existência e unicidade do P.V.I. Desta

forma basta demonstrar a função y definida acima é a solução. Pela condição inicial, vemos que
y(0) = Eα (0) = 1, visto que
λxα λ2 x2α
Eα (λxα ) = 1 + + + ··· ;
Γ (1 + α) Γ (1 + 2α)
além de que, no caso de n ≥ 2, isto é, α > 1, temos y (k) = 0, para k = 1, 2, . . . , n − 1, desde que
λxα λ2 x2α
y(x) = 1 + + + ··· ;
Γ (1 + α) Γ (1 + 2α)
o que implica que
λxα−k λ2 x2α−k
y (k) (x) = + + ··· ;
Γ (1 + α − k) Γ (1 + 2α − k)
Com relação ao operador diferencial, vamos considerar o caso λ = 0, segue que y(x) = Eα (0) =
1 e consequentemente C Daα y(x) = 0 = λy(x). Se λ ̸= 0 segue que
" ∞
#
α
C
X (λx )
Daα y(x) = C α
Da (x)
j=0
Γ (1 + jα)
"∞ #
X (λj xαj )
= I0n−α Dn (x)
j=0
Γ (1 + jα)
"∞ #
X (λj Dn xαj )
= I0n−α (x)
j=0
Γ (1 + jα)
"∞ #
X (λj Dn xαj )
= I0n−α (x)
j=1
Γ (1 + jα)
"∞ #
X (λj xαj−n )
= I0n−α (x)
j=1
Γ (1 + jα − n)
∞
X (λj I n−α xαj−n )
0
= (x)
j=1
Γ (1 + jα − n)
∞
X (λj I n−α xαj−n )
0
= (x)
j=1
Γ (1 + jα − n)
∞
X (λj xαj−α )
= (x)
j=1
Γ (1 + jα − α)
∞
X (λj+1 xαj )
= (x)
j=0
Γ (1 + jα)
∞
X (λxα )j
= λ (x) = λy(x)
j=0
Γ (1 + jα)

63
5.2 Métodos de Diferenças Fracionárias Regressivas

onde usamos as propriedades de convergências da série que define a função de Mittag-Leffler

5.2 Métodos de Diferenças Fracionárias Regressivas

Nesta seção iremos desenvolver o algoritmo para a solução da equação diferencial fracionais em
três contextos. Primeiro iremos abordar o operador diferencial tipo Grünwald-Letnikov,
GL α
Daα y(x) = f (x, y(x)), Dα−k y(0) = 0, 0 < α < 1 (5.13)
Seguimos com a equação do tipo Riemann-Liouville.
Dα y(x) = f (x, y(x)), Dα−k y(0) = bk , k = 1, 2, . . . , n − 1, lim I n−α y(z) = bn (5.14)
z→0+
E por último, abordamos o operador diferencial de Caputo.
C
Daα y(x) = f (x, y(x)), Dk y(0) = bk , k = 0, 2, . . . , n − 1, (5.15)
onde α > 0, α ∈ / N e n = ⌈α⌉. Nosso interesse é a solução de y(x) num intervalo fechado [0, X].
Seja a malha de pontos x0 , x1 , . . . , xN , onde N = X/h, x0 = 0 e xk = x0 + kh. Denotamos ym
a aproximação de y(xm ) e igualmente fm = f (xm , ym ). Iremos desenvolver a versão não inteira
do conhecido método clássico de diferenças regressivas. Existem diversas formas de operadores
diferenciais não inteira. Iremos iniciar com o operador diferencial de Grünwald-Letnikov.

5.2.1 Diferenças Regressivas e operador Grünwald-Letnikov

Na equação (5.13) temos as condições iniciais homogênea e como resultado do Corolário 4.4,
temos que o problema (5.13) é equivalente aos problemas (5.14) e (5.15) para 0 < α < 1. A vantagem
da formulação (5.13) é que obtemos uma técnica de discretização imediata a partir da definição do
Derivada de Grünwald-Letnikov.

l
(∆αh y)(x)

GL 1 X k α
Daα y(x) = lim = lim α (−1) y(x − kh), lh = x, α > 0
h→0 hα h→0 h
k=0
k
Numericamente, adotamos a forma discreta do operador, para um valor h ≪ 1.
m
GL 1 X k α
Daα F y(xm ) = α (−1) y(xm − kh), m = 0, 1, . . . , N,
h k=0 k

na malha de pontos x0 , x1 , . . . , xN . Desta forma a discretização o problema (5.13).

m
1 X
ωk y(xm − kh) = f (xm , y(xm )), m = 0, 1, . . . , N,
hα k=0

onde
αk
ωk = (−1) .
k

64
5.2 Métodos de Diferenças Fracionárias Regressivas

Podemos resolver este sistema de equações em cada ponto xm da malha por

X m
α
ym = h f (xm , ym ) − ωk y(xm − kh), m = 1, 2, . . . , N y0 = 0 (5.16)
k=1
Este conjunto de equações fornece a solução numérica de (5.13). No cálculo de cada passo m, ym
aparece de forma implı́cita, mas sendo a única variável desconhecida, pois y1 , y2 , . . . ym−1 foram
obtidos em passos anteriores.
Em geral, teremos um sistema de equações não lineares e portanto teremos que aplicar o método
de ponto fixo.
A versão finita do operador Grünwald-Letnikov fornece uma aproximação de primeira ordem do
operador de diferencial de Riemann-Liouville e pela escolha da condição inicial, também aproxima o
operador de Caputo com ordem um. Assim, a fórmula (5.13) nos dá um método numérico de primeira
ordem para resolver equações do tipo (5.13), assim como as equações do tipo (5.14) e (5.15), quando
0 < α < 1 e condição inicial y(0) = 0. Antes de generalizar o problema (5.13), para os casos α > 1
e condições iniciais não homogênea, vamos analisar o pesos ωk .
Os coeficientes ωk pode ser calculados de forma recursiva, sendo ω0 = 1, para k ∈ N, segundo
a Definição 3.8 segue que

k α
ωk = (−1) (5.17)
k
(−1)k−1 Γ (k − α)
= (−1)k
Γ (k + 1)Γ (1 − α)
−(−1)k−2 Γ (k − 1 − α)(k − 1 − α)
= −(−1)k−1
kΓ (k)Γ (1 − α)
k−2
(−1) Γ (k − 1 − α) k − (α + 1)
= (−1)k−1
Γ (k)Γ (1 − α) k

α + 1)
= 1− ωk−1
k
(5.18)
para todo k ∈ N. Outra forma de se calcular ωk é pela generalização da função
ω(ζ) = (1 − ζ)α , (5.19)
isto é, os k primeiros coeficientes da Série de Taylor de (5.19) são os k primeiros pesos ωk .
A princı́pio, temos que a fórmula (5.16) converge com ordem um. Este fato é mal interpretado
que (5.16) comporta-se como O(h) quando h → ∞. Precisamos de um peso adicional inicial em
cada etapa, dada por
m−α Γ (α)
ωm = − (−1)m (5.20)
Γ (1 − α) Γ (α − m)Γ (m + 1)
m
m−α X
= − ωj , )
Γ (1 − α) j=0

65
5.2 Métodos de Diferenças Fracionárias Regressivas

resulta na formulação modificada

m m
!
α
X m−α X
ym = h f (xm , ym ) − ωk y(xm − kh) − − ωj , y0 , m = 1, 2, . . . , N, (5.21)
k=1
Γ (1 − α) j=0
cujo o erro é da ordem O(h). Observe que as formulações (5.16) e (5.21) são equivalentes no caso de
condições iniciais homogênea.
Vamos considerar a classe de problemas que vamos aplicar este esquema de diferenças para o
problema (5.13) com α > 0 e correspondente condições iniciais na forma de Riemann-Liouville ou
Caputo. Este tipo de problemas são, em parte, abordado no trabalho de Podlubny [35], usando uma
abordagem matricial.
A ideia básica da abordagem de Podlubny, é baseada no fato que a equação diferencial com
condições iniciais homógenas (5.16) pode ser representado na forma matricial
    
ω0 y0 f0
ω ω0   y 1   f1 
    
1 
 .1 α
  ..  =  ..  ⇔ BN YN = FN (5.22)
..     
hα  . ..
 . . .  .   . 
ωN ωN −1 · · · ω0 yN fN
Sendo a condição inicial homogênea, isto é y0 , a primeira equação pode ser descartada. Obtemos a
solução nos pontos da malha resolvendo o sistema linear, em vez de usar a forma recursiva (5.16).
No entanto essa abordagem apresenta algumas restrições. Podlubny descreve o método para equações
diferenciais não inteira de k termos, ou seja, equações com mais de um operador diferencial com
condições iniciais homogêneas ou não homogêneas do tipo Riemann-Liouville ou Caputo.
Descrevendo o método para equações de k termos, os tipos de problemas que podem ser resolvidos
é ampliado e a ideia básica é discretizar os operadores diferenciais pela fórmula (5.16). Por outro
lado, a restrição a problemas lineares é necessária e um limitante para aplicação do método.
No problema do tipo Riemann-Liouville, as condições iniciais não são do tipo clássico, ou seja,
usa informações de derivada não inteira no ponto inicial x0 , que muitas vezes não tem significado
fı́sico conhecido ou pode ser medido. No caso do tipo Caputo, as condições iniciais são de ordem
inteira. Assim sendo, vamos analisar soluções numéricas para o problema (5.15) do tipo Caputo.
Segundo o Corolário 4.4, o operador Grünwald-Letnikov e o operador de Caputo se relacionam
da seguinte forma
C α
Da y(x) = GL Daα y(x) − Dα Tn−1 [y; 0](x)

onde y(x) é assumido ser n-vezes continuamente diferenciáveis. Aplicado ao problema (5.15), o
polinômio de Taylor Tn−1 [y; 0](x), é completamente definido pelas condições iniciais e segue que, a
equação (5.21) no caso α > 0 pode ser reescrita como
m
!
α
X m−α
ym = h f (xm , ym )− ωk y(xm −kh)− y0 +hα Dα Tn−1 [y; 0](xm ), m = 1, 2, . . . , N
Γ (m − α) − m
P
k=1 j=0 ωj
(5.23)
Com isso obtemos um método de primeira ordem, do tipo Caputo, que soluciona equações lineares
ou não lineares, homogêneas ou não homogêneas para α > 0.

Observação

66
5.2 Métodos de Diferenças Fracionárias Regressivas

a) Observamos que a derivada Dα Tn−1 [y; 0](x) na fórmula (5.23) pode ser calculada por
n−1
α
X bk xk−α
m
D Tn−1 [y; 0](xm ) =
k=0
Γ (k + 1 − α)
b) Ao fazer a análise da equação (5.23) para os casos não lineares de Caputo, observamos que o
polinômio de Taylor é nulo para as condições iniciais homogêneas. Logo o sistema de equações
(5.22) permanece o mesmo para α > 1 e 0 < α ≤ 1.
c) Ao considerar equações do tipo Riemann-Liouville (5.14), obtemos uma fórmula semelhante a
(5.23). Isto mostra a conexão entre os operadores de Riemann-Liouville e Grünwald-Letnikov.
Podemos usar a seguinte identidade
Dα y(x) = Dα (y(x) − T̂n−1 [y; 0](x)),
onde n
X bk xα−k
T̂n−1 [y; 0](x) =
k=1
Γ (α − k + 1)

é a expansão de Taylor no contexto do operador Riemann-Liouville.

Resumindo o método numérico baseado na diferença finita fracionária:
1. Uma equação diferencial fracionária linear ou não linear, do tipo Riemann-Liouville ou Caputo,
com α > 0 e condições homogêneas podem ser resolvido por (5.16), onde, em cada etapa
m = 1, . . . , N , uma equação linear ou não linear pode ser resolvida.
2. Uma equação diferencial fracionária linear, do tipo Riemann-Liouville ou Caputo com α > 0
com condições iniciais não homogêneas pode ser convertido num problema com condições
homogêneas pela transformação
Xn−1 Xn
k
y(x) = bk x + z(x) ou y(x) = bk xα−k + z(x)
k=0 k=0
para o tipo Caputo ou Riemann-Liouville.
Para uma equação diferencial fracionária não linear do tipo Caputo, com α > 0 e condições iniciais
não homogêneas pode ser resolvida por (5.21).
3. Para uma equação diferencial fracionária não linear do tipo Riemann-Liouville, com α > 0 e
condições iniciais não homogêneas pode ser resolvida por nenhum dos métodos apresentados.
(5.21).

5.2.2 Método de diferenças finitas regresivas de Lubich

Nesta seção será apresentado um método de diferenças finitas de alta ordem baseado no método
clássico. A ideia foi apresentada por Lubich [24–27].
Vamos considerar a equação integral de Abel-Volterra
Z x
1
y(x) = f (x) + (x − t)α−1 K(x, t, y(t))dt, x ∈ [0, X], α > 0 (5.24)
Γ (α) 0
com o kernel K limitado e com termo forçante f (x). O método multipassos linear fracionário é dado
por:

67
5.2 Métodos de Diferenças Fracionárias Regressivas

m
X X
α
ym = f (xm ) + h ωm−j K(Xm , xj , y(xj )) + hα j = 0m wm,j K(Xm , xj , y(xj )) (5.25)
j=0

que fornece uma aproximação para y(xm ), com erro satisfazendo

max |ym − y(xm )| = O(hp−ϵ ), (5.26)
0≤m≤N
onde ϵ > 0 e pequeno. Os coeficientes da convolução, em 5.25, ωm são gerados pela função
α
α σ(1/ζ)
ω = ,
ρ(1/ζ)
onde (ρ, σ) são os polinômios caracterı́stico do método multipasso linear clássico e os coeficientes
wm,j são dados pelo sistema linear
s m
X
γ Γ (1 + σ) α+σ
X
wm,j j = m − ωm−j j γ , γ ∈ A (5.27)
j=0
Γ (1 + σ + α) j=1
com
A = {γ = k + jα; k, j ∈ N0 , γ ≤ p − 1}, #A = s + 1 (5.28)

A equação diferencial fracionária do tipo Caputo (5.15) pode ser interpretada como uma equação
integral de Abel-Volterra
Z x
1
y(x) = Tn−1 [y; 0](xm ) + (x − t)α−1 f (t, y(t))dt (5.29)
Γ (α) 0
onde o termo forçante Tn−1 [y; 0](xm ) é completamente definido pela condição inicial de (5.15) e o
kernel K(t, x, y(x)) = f (x, y(x)) é dado pelo termo do lado direito de (5.15). O método diferenças
regressivo clássico de ordem p possui função geradora
p p
X
k
X 1
ω(ζ) = ωk ζ = (5.30)
k=0 k=1
k(1 − ζ)k
Com estes resultados, podemos descrever o método de diferenças regressivas fracionária.
Teorema 5.5
Seja α > 0 e n = ⌈α⌉. O método de diferenças regressivo fracionário de Lubich, de ordem
p ∈ {1, 2, . . . , 6}, para uma equação do tipo Caputo (5.15) reescrita como uma equação de
Abel-Volterra, é dada por
m
X s
X
α α
ym = Tn−1 y; 0(xm ) + h ωm−j f (xj , y(xj )) + h wm,j f (xj , y(xj )) (5.31)
j=0 j=0
para m = 1, . . . , N , onde os coeficientes ωm são obtidos pela função geradora
p
!−α
X 1
ω α (ζ) = (1 − ζ)k (5.32)
k=1
k
e com os coeficientes iniciais wm,j são obtidos na solução do sistema linear (5.27). A equação
(5.31) fornece uma aproximação de ordem O(hp−ϵ ) com ϵ > 0 para todos os pontos da malha
xm . ♡

No próximo teorema damos uma versão modificada do método.

68
5.3 Método do tipo Adams-Bashforth-Moulton

Teorema 5.6
Seja α > 0 e n = ⌈α⌉. O método de diferenças regressivo fracionário de Lubich, de ordem
p ∈ {1, 2, . . . , 6}, para uma equação do tipo Caputo (5.15) é dada por é dada por
m−1
X Xs
α
ym = h f (xm , ym ) − ωm−j y(xj ) − wm,j y(xj ) + hα Dα Tn−1 [y; 0](xm ) (5.33)
j=0 j=0
para m = 1, . . . , N , onde os coeficientes wm,j são obtidos pela função geradora
s m
X
γ Γ (1 + σ) α+α
X
wm,j j = m − ωm−j j γ , γ ∈ A (5.34)
j=0
Γ (1 + σ + α) j=1

com A definido em (5.28). O coeficientes wm,j tem ordem O(mα−1 ) e o erro satisfaz O(hp−ϵ )
com ϵ > 0 para todos os pontos da malha xm . ♡

Observação
a) A restrição de que a ordem p ∈ {1, 2, . . . , 6} se deve ao fato de que o método clássico é estável
até a ordem p = 6.
b) O valor de ϵ no Teorema 5.5 e no Teorema 5.6 é dado pela equação
ϵ = p − α − min{γ = k + jα, k, j ∈ N0 , γ > p − 1}
c) A vantagem da formulação (5.34) comparada com a formulação (5.32) é que na formulação (5.34)
os pesos são calculados em função dos pesos já obtidos.
d) O Teorema 5.5 e o Teorema 5.6 podem ser formulados para equações do tipo Riemann-Liouville.
e) O Teorema 5.5 e o Teorema 5.6 iniciam os coeficientes de forma particular. Garantem que o
conjunto de funções base b(x) = xk+jα , k + jα ∈ A são integráveis ou deriváveis de forma
exata.
O Teorema 5.6 estabelece uma generalização do método de diferenças finitas clássico para o caso
fracionário. Ainda não abordamos os aspectos computacionais do cálculo dos coeficientes ωm e Wm,j
e o fato da formulação (5.34) não pode ser resolvidos passo a passo para os primeiros s + 1 pontos da
malha. Todavia o cálculo dos s + 1 iniciais, pode ser feito de forma simultânea, com a resolução de
um sistema nâo linear.

5.3 Método do tipo Adams-Bashforth-Moulton

Nesta seção iremos apresentar uma versão do método de Adams-Bashforth-Moulton para resolver
uma equação diferencial fracionária Em particular iremos considerar o problema na forma da integral
Abel-Volterra.
n−1 k Z x
X x 1
y(x) = bk + (x − t)α−1 f (t, y(t))dt (5.35)
k=0
k! Γ (α) 0

Suponha que (5.35) tenha solução única para x ∈ [0, X]. Estamos interessados em uma solução
numérica na malha uniforme M = {xj = jh : j = 0, 1, 2, . . . , N } para algum N inteiro positivo e
h = X/N .
Uma caracterı́stica dos métodos de Adams-Bashfort é que este é um método de passo múltiplo, ou

69
5.3 Método do tipo Adams-Bashforth-Moulton

seja, necessita de um método de passo simples de mesma ordem para obter as primeiras aproximações
yj ≈ y(xj ), j = 1, 2, 3, . . . , k. Além disso, podem ser aplicados no esquema implı́cito, ou seja, o
cálculo de yk+1 depende do conhecimento de uma aproximação ỹk+1 também obtido por um método
de passo simples, configurando um esquema Preditor-Corretor.
Inicialmente vamos considerar que temos as primeira aproximações calculadas. Usaremos uma
formulação da quadratura trapezoidal para aproximar a integral de (5.35), onde os pontos da mala
xj , j = 1, 2, . . . , k + 1 são tomados com respeito a função peso (xk+1 − ·)α−1 . Em outras palavras,
aplicamos a aproximação
Z xk+1 Z xk+1
α−1
(xk+1 − z) g(z)dz ≈ (xk+1 − z)α−1 g̃k+1 (z)dz (5.36)
0 0
onde g̃k+1 é a interpolação linear de g nos pontos da malha xj , j = 1, 2, . . . , k + 1.
Usaremos um esquema do tipo Newton-Côtes para aproximar a integral (5.36).

Z xk+1 k+1
X
α−1
(xk+1 − z) g̃k+1 (z)dz = aj,k+1 g(xj ) (5.37)
0 j=0
Z xk+1
aj,k+1 = (xk+1 − z)α−1 ϕj,k+1 (z)dz (5.38)
0

 (z − xj−1 )/(xj − xj−1 )
 se xj−1 < z ≤ xj ,
ϕj,k+1 (z) = (xj+1 − z)/(xj+1 − xj ) se xj < z ≤ xj+1 , (5.39)

0

caso contrário
Com o cálculo da integral (5.38) para um xj obtemos
(xk+1 − x1 )α+1 + xαk+1 (αx1 + x1 − xk+1
a0,k+1 = (5.40)
x1 α(α + 1)
(xk+1 − xj−1 )α+1 + (xk+1 − xj )α (α(xj−1 − xj ) + xj−1 − xk+1 )
aj,k+1 =
(xj − xj−1 )α(α + 1)
α+1
(xk+1 − xj+1 ) − (xk+1 − xj )α (α(xj − xj+1 ) − xj+1 + xk+1 )
+ se 1 ≤ j ≤ k
(5.41)
(xj+1 − xj )α(α + 1)
(xk+1 − xk )α
ak+1,k+1 = (5.42)
α(α + 1)
No caso de uma malha equidistante, (xj = jh para um h fixo) então estas relações se reduzem a

hα

k α+1 − (k − α)(k + 1)α


 se j = 0,



 α(α + 1)



hα


(k − j + 2)α+1 + (k − j)α+1 − 2(k − j + 1)α+1

aj,k+1 = se 1 ≤ j ≤ k, (5.43)
 α(α + 1)





hα




 se j = k + 1.
α(α + 1)

Com isto temos a variante fracional para o esquema corretor de um passo dado por:
n−1 j k
!
X xk+1 1 X
yk+1 = bj + aj,k+1 f (xj , yj ) + ak+1,k+1 f (xk+1,yk+1
P ) (5.44)
j=0
j! Γ (α) j=0

70
5.4 Exemplos Numéricos

onde ykP é o valor obtido pelo esquema preditor, que iremos abordar a seguir.
A ideia é usar a generalização do esquema de Adams-Bashforth de um passo. usando a mesma
técnica aplicada no caso corretor. Usamos a regra de quadratura para co cálculo da integral (5.35)
Z xk+1 K
X
(xk+1 − z)α−1 g(z)dz ≈ bj,k+1 g(xj ), (5.45)
0 j=0
onde
α)
xj+1 (xk+1 − xj )α − (xk+1 − xj+1
Z
α−1
bj,k+1 = (xk+1 − z) dz = . (5.46)
xj α
Estamos considerando aproximação constantes por partes e não linear por partes, e usamos a
Função de Haar em [xj , xj+1 ]. Novamente, no caso equidistante, temos a expressão mais simples
hα
bj,k+1 = ((k + 1 − j)α − (k − j)α ) (5.47)
α
P
Desta forma, o esquema preditor yk+1 é determinado pelo esquema Adams-Bashforth
n−1 j k
P
X xk+1 1 X
yk+1 = bj + bj,k+1 f (xj , yj ) (5.48)
j=0
j! Γ (α) j=0
Em resumo, o esquema Adams-Bashforth-Moulton fracionário é completamente descrito pelas equações
(5.48) e (5.44), onde os coeficientes aj,k+1 e bj,k+1 são definidos em (5.46) e (5.43).

5.4 Exemplos Numéricos

Seja o Problema de Valor Inicial dado por
(
C α 2
Da y = t2 + Γ (3−α)
t2−α −y
(5.49)
y(0) = 0
cuja a solução exata é dada por y(t) = t2 , para 0 < α ≤ 1. Usaremos o esquema numérico de esquema
Adams-Bashforth-Moulton fracionário, (Ver Apêndice C), definido pelas equações equações (5.48) e
(5.44), onde os coeficientes aj,k+1 e bj,k+1 são definidos em (5.46) e (5.43).
Testamos o esquema, usando uma malha de 1000 pontos para aproximar a solução da EDOF no
intervalo de tempo [0, 1] e α = 0.5. A aproximação e o erro absoluto são apresentados na Figura 5.1.
Um outro exemplo é apresentado para o Problema de Valor Inicial dado por
(
Da y = Γ (5+α)
C α
24
∗ t4 − y 2 + t8+2∗α
(5.50)
y(0) = 0
cuja a solução exata é y(t) = t4 para 0 < α ≤ 2.
Na Figura 5.2 mostra a aproximação e o erro onde foi usado N = 1000, α = 1.5 e t ∈ [0, 1].
Neste caso, em que α = 1.5, o número de condições iniciais é dado por ⌈α⌉ = 2. No exemplo
apresentado, esses valores são nulos.
Para exemplificar a aplicação, num sistema de equações diferenciais fracionárias acopladas, usa-
remos a Equação de Bagley-Torvik [49], onde a derivada fracionária surge naturalmente na descrição
de certos movimentos de um fluido newtoniano.

71
5.4 Exemplos Numéricos

(a) Solução exata e aproximada (b) Erro absoluto

Figura 5.1: Aproximação com N = 1000, no intervalo de tempo [0, 1] e α = 0.5

(a) Solução exata e aproximada (b) Erro absoluto

Figura 5.2: Aproximação com N = 1000, no intervalo de tempo [0, 1] e α = 1.5

3/2
Ay ′′ (t) + B C D0 y(t) + Cy(t) = f (t) (5.51)

onde A, B e C são constantes reais. A equação (5.51) pode ser reescrita como um sistema de
quatro equações diferenciais acopladas de ordem α = 1/2.
 C 1/2

 D0 y1 (t) = y2 (t)
C 1/2

 D0 y2 (t) = y3 (t)
C 1/2
(5.52)


 D0 y3 (t) = y4 (t)
C 1/2
A−1 (−Cy1 (t) − By4 (t) + f (t))

D0 y4 (t) =
′
com condições iniciais y1 (0) = y0 , y2 (0) = 0, y3 (0) = y0 , y4 (0) = 0
Vamos considerar o caso em que A = B = C = 1 e f (t) é dada por
t3/2
f (t) = t3 + 7t + 1 + 8 √
π
Neste caso a equação diferencial apresenta a solução exata y(t) = t3 + t + 1.
Testamos o esquema, usando N = 1000 para aproximar a solução da EDOF no intervalo de

72
5.4 Exemplos Numéricos

tempo [0, 6] e α = 0.5. A aproximação e o erro absoluto são apresentados na Figura 5.3.

(a) Solução exata e aproximada (b) Erro absoluto

Figura 5.3: Aproximação com N = 1000, no intervalo de tempo t ∈ [0, 6] e α = 0.5

Além da solução da equação diferencial, o procedimento acima também fornece aproximações

das derivadas fracionária de ordem α = 1/2, 1, 3/2 da solução. Na Figura 5.4 apresentamos estas
aproximações.

Figura 5.4: Derivadas de ordem α = 0, 1/2, 1, 3/2, 2 com N = 1000, no intervalo de tempo [0, 6]

Vamos apresentar um outro exemplo, onde consideramos uma equação diferencial, onde a ordem
é inteira e vamos transformá-la num sistema de equações diferenciais fracionárias e ver como se
comporta os procedimentos numéricos. Para isto vamos considerar o P.V.I. de segunda ordem

′′
 y +y =0

y(0) = 0 (5.53)

 ′
y (0) = 1
cuja a solução exata é dada por y(t) = sen(t). Transformando o P.V.I. num sistema de quatro equações

73
5.4 Exemplos Numéricos

diferenciais fracionárias de ordem α = 1/2 temos que

 C 1/2

 D0 y1 (t) = y2 (t)
 C D1/2 y (t) =

y3 (t)
0 2
C 1/2 (5.54)


 D0 y3 (t) = y4 (t)
 C 1/2
D0 y4 (t) = −y1 (t)
com condições iniciais y1 (0) = 0, y2 (0) = 0, y3 (0) = 1, y4 (0) = 0
Testamos o esquema, com N = 1000, para t ∈ [0, 6]. Na Figura temos a comparação da solução
exata e a aproximação e o erro 5.5 absoluto cometido.

(a) Solução exata e aproximada (b) Erro absoluto

Figura 5.5: Aproximação com N = 1000, no intervalo de tempo t ∈ [0, 6] e α = 0.5

Como no exemplo da Equação de Bagley-Torvik, podemos ver como se comporta as derivadas

fracionárias da solução, representadas na Figura 5.6.

74
5.4 Exemplos Numéricos

Figura 5.6: Derivadas de ordem α = 0, 1/2, 1, 3/2, 2 com N = 1000, no intervalo de tempo [0, 1]

75
Capı́tulo 6 Equação Diferencial Parcial
Fracionáriado tipo Riesz
Vários modelos nas áreas de finanças, engenharia e fı́sica, tem feito uso das equações diferenciais
parcial fracionária. As equações diferenciais parciais fracionárias com derivada de Riesz têm sido
amplamente utilizadas para modelar fenômenos de difusão anômala e transporte em meios complexos.
Nesta capı́tulo iremos abordar a formulação matemática dessas equações, com ênfase nas propriedades
da derivada fracionária de Riesz [52]. Apresentamos três esquemas numéricos no contexto da derivada
fracionária de Riesz. Vamos considerar a equação cinética fracionária, com condições de fronteiras
de Dirichlet homogêneas.

6.1 Equação de Advecção-Difusão Fracionária do tipo de Riesz

Como modelos de uma equação diferencial parcial fracionária, vamos considerar a equação de
advecção-difusão fracionária do tipo de Riesz, com condições de fronteiras de Dirichlet homogêneas
dada pela definição:
Definição 6.1
Sejam 1 < α ≤ 2, 0 < β ≤ 1, 0 < t < T e 0 < x < L. Então
∂u(x, t) ∂α ∂β
= Kα u(x, t) + K β u(x, t) (6.1)
∂t ∂|x|α ∂|x|β
sujeita as condições
u(0, t) = u(L, t) = 0 (6.2)
u(x, 0) = g(x), (6.3)
é chamada de Equação de Advecção-Difusão Fracionária do tipo Riesz (EADFR ), onde
derivadas fracionárias do espaço de Riesz são definidas por:
∂α α α
α
u(x, t) = −Cα (D[a,x] + D[x,a] ), (6.4)
∂|x|
∂β β β
β
u(x, t) = −Cβ (D[a,x] + D[x,a] ), (6.5)
∂|x|
com
1 1
Cα = , Cβ =
2 cos(πα/2) 2 cos(πβ/2)
♣
α α
As definições de D[a,x] e D[x,a] são dadas na Observação 4.1. Como exemplo, u pode representar
a concentração de soluto onde Kα e Kβ representam o coeficiente de dispersão e velocidade média
do fluı́do respectivamente. No contexto fı́sico do modelo de transporte advecção-difusão fracionária
vamos assumir que Kα > 0 e Kβ ≥ 0 tal que o fluxo é da esquerda para a direita.
6.1 Equação de Advecção-Difusão Fracionária do tipo de Riesz

6.1.1 Laplaciano Fracionário

Para verificar a existência de soluções analı́ticas e numéricas, vamos precisar de alguns resultados
preliminares sobre o Laplaciano Fracionário. Iniciamos com a sua definição.
Definição 6.2 (Laplaciano Fracionário)
Suponha que o Laplaciano (−∆) tenha um conjunto completo ortogonal de autofunções φn
correspondentes aos autovalores λ2n numa região limitada D, isto é,(−∆)φn = λ2n φn , na região
D; B(φ) = 0 em ∂D, onde B(φ) representa uma das três condições homogêneas de fronteira.
Seja
( +∞ +∞
)
X X
G= f= cn φn , cn = ⟨f, φn ⟩, |cn |2 |λn |γ < ∞, γ = max(α, 0) .
n=1 n=1
Então para qualquer f ∈ G, definimos o Laplaciano Fracionário por
+∞
X
α/2
−(−∆) f = cn (λ2n )α/2 φn
n=1
♣

O lema que segue mostra a relação entre o Laplaciano Fracionário e a derivada fracionária no
espaço de Riesz.
Lema 6.1
Para uma função y(x) definida em R que tende a zero quando x → ±∞, vale:
∂α
−(−∆)α/2 y(x) = y(x)
∂|x|α
onde −(−∆)α/2 é o operador Laplaciano fracionário definido via Transformada de Fourier.
♡

Demonstração Segundo Yang [Link] [52], o operador Laplaciano fracionário é definido como
−(−∆)α/2 y(x) = −F −1 |x|α Fy(x),
onde F e F −1 representam a Transformada de Fourier e sua inversa. Consequentemente temos
Z +∞ Z +∞
α/2 1 −iξx α iξη
−(−∆) y(x) = − e |ξ| e u(η)dη dξ
2π −∞ −∞
Supondo que y(x) se anula em x → ±∞, aplicamos a integração por partes
Z +∞
1 +∞ iξη ′
Z
iξη
e y(η)dη = e y (η)dη
−∞ iξ −∞
Com isto obtemos
Z +∞ Z +∞ α

α/2 1 ′ i(η−ξ) |ξ|
−(−∆) y(x) = − y (η) i e dξ dη
2π −∞ −∞ ξ
Vamos chamar de I o termo entre parênteses da equação acima. Então
Z +∞ α
Z +∞ Z +∞
i(η−ξ) |ξ| i(ξ−η) α−1 i(η−ξ) α−1
I=i e dξ = i − e ξ dξ + e ξ dξ
−∞ ξ 0 0
Note que, sendo
Z +∞
Γ (v)
p−1
L(t ) = e−st tp−1 dt = p , Re(v) > 0
0 s

77
6.1 Equação de Advecção-Difusão Fracionária do tipo de Riesz

e 0 < α < 1, temos que

−Γ (α) Γ (α) sign(ξ − η)Γ (α)Γ (1 − α) α−1 α−1

I=i + = t + (−i)
i(η − ξ) i(ξ − η) |ξ − η|α Γ (1 − α)
π
Usando os fatos Γ (α)Γ (1 − α) = e iα−1 + (−i)α−1 = 2sen(πx/2), obtemos
sen(πx)
sign(ξ − η)π
I= .
cos(πx/2)|ξ − η|α Γ (1 − α)
Por isso, para 0 < α < 1, segue

+∞
sign(ξ − η)π
Z
1
−(−∆) α/2
y(x) = − u′ (η) dη
2π
−∞ cos(πx/2)|ξ − η|α Γ (1 − α)
Z +∞ Z +∞
u′ (η) u′ (η)

1 1 1
= − dη − dη .
2 cos(πx/2) Γ (1 − α) −∞ (ξ − η)α Γ (1 − α) x (ξ − η)α
Segundo o Teorema 4.20, para 0 < α < 1, a derivada fracionária de Grünwald-Letnikov em [a, x] é
dada por Z x
GL α y(a)(x − a)−α 1 y ′ (η)
D [y(x)][a,x] = + dη
Γ (1 − α) Γ (1 − α) a (x − η)α
Desta forma, se y(x) tende a zero, quando a → −∞, temos que
x
y ′ (η)
Z
GL α 1
D [y(x)][−∞,x] = dη
Γ (1 − α) −∞ (x − η)α
De forma análoga, se y(x) tende a zero, quando b → +∞, temos que
+∞
y ′ (η)
Z
GL α 1
D [y(x)][x,+∞] =− dη.
Γ (1 − α) x (η − x)α
Se y(x) é contı́nua e y ′ (x) é integrável para x ≥ a, então para todo α no intervalo (0, 1) a derivada
de Riemann-Liouville existe e coincide com a a derivada de Grünwald-Letnikov. Finalmente, para
0 < α < 1 temos
1 ∂α
−(−∆)α/2 y(x) = −
GL α
D [y(x)][−∞,x] + GL Dα [y(x)][x,+∞] ) =

y(x)
2 cos(πα/2) ∂|x|α
Para o caso em que 1 < α < 2, a demostração pode ser feita de forma similar a apresentada para
o caso 0 < α < 1. [52]
E com isto, podemos concluir que sendo n − 1 < α < n então

1 ∂α
−(−∆)α/2 y(x) = −
GL α
D [y(x)][−∞,x] + GL Dα [y(x)][x,+∞] ) =

y(x),
2 cos(πα/2) ∂|x|α
onde Z x
GL α 1 ∂n y(η)
D [y(x)][−∞,x] = dη,
Γ (n − α) ∂x −∞ (x − η)α+1−n
n
Z +∞
GL (−1)n ∂ n y(η)
Dα [y(x)][x,+∞] = n
dη
Γ (n − α) ∂x x (η − x)α+1−n

No cado de funções y(x) definidas num intervalo finito [0, L], os resultados acima são válidos

78
6.2 Esquemas Numéricos

para a função estendida

(
y(x), x ∈ (0, L)
ŷ(x) =
0, x ∈ (−∞, 0] ∪ [L, +∞)

6.1.2 Solução Analı́tica

A representação espectral do operador Laplaciano (−∆) é obtida por intermédio da resolução do
problema de autovalor-autofunção
−∆φ = λφ
φ(0) = φ(L) = 0.
π 2 n2
cuja a solução fornece os autovalores λn = 2 para n = 1, 2, . . . , e as correspondentes autofunções
nπ L
são φn = sen .
L
Vamos considerar que a solução de EADFR definido em (6.1) seja da forma
X∞
u(x, t) = cn (t)φn (x)
n=1
que satisfaz as condições de fronteiras, independente de cn (t). Usando a Definição 6.2 e substituindo
u(x, t) na equação diferencial (6.1), obtemos
∞
X d α/2 β/2

cn (t) + Kα (λn ) + Kβ (λn ) cn (t) φn (x) = 0
n=1
dt
Com isto segue que os coeficientes cn (t) devem satisfazer a equação diferencial ordinária
d
cn (t) + Kα (λn )α/2 + Kβ (λn )β/2 cn (t) = 0

dt
cuja a solução geral é da forma
c (t) = c (0)e−[Kα (λn ) +Kβ (λn ) ]t
α/2 β/2
n n

Para determinar cn (0) usamos as condição inicial definida em (6.3)

+∞
X
u(x, 0) = cn (0)φn (x) = g(x)
n=1
da qual se deduz Z L
2 nπη
cn (0) = g(η)sen dη = bn
L 0 L
Desta forma, temos a solução analı́tica de (6.1) dada por
+∞ nπ α/2 β/2
e−[Kα (λn ) +Kβ (λn ) ]t
X
u(x, t) = bn sen (6.6)
n=1
L

6.2 Esquemas Numéricos

Nesta seção iremos apresentar três esquemas numéricos para EADFR descrita na Definição 6.1.

79
6.2 Esquemas Numéricos

6.2.1 Esquema de aproximação L1-L2

Nesta seção iremos abordar os esquemas numéricos baseados na Definição 6.1 e no Lema 6.1.1.
Vamos iniciar com a discretização do termo de difusão. Segundo o Teorema 4.20, para 1 < α ≤ 2 no
intervalo finito [0, L] temos os operadores diferenciais Grünwald-Letnikov a esquerda e direita dadas
por:

u(0)x−α u′ (0)x1−α x
u(2) (η)
Z
GL α 1
D [u(x)][0,x] = + + α−1
dη, (6.7)
Γ (1 − α) Γ (2 − α) Γ (2 − α) 0 (x − η)
u(L)(L − x)−α u′ (L)(L − x)1−α
Z L
GL 1 u(2) (η)
Dα [u(x)][x,L] = + + dη. (6.8)
Γ (1 − α) Γ (2 − α) Γ (2 − α) x (η − x)α−1

A relação entre as definições de Riemann–Liouville e Grünwald–Letnikov (Ver Teo. 4.21 )

permite o uso da definição de Riemann–Liouville durante a formulação do problema e, em seguida,
da definição de Grünwald–Letnikov para a obtenção da solução numérica. Portanto, usando a ligação
entre essas duas definições e o algoritmo proposto por [Oldhann74], podemos facilmente obter o
seguinte esquema de discretização numérica para o termo de difusão.
Suponha que o domı́nio espacial é [0, L]. A malha composta por N intervalos de h = L/N
e xl = lh, para 0 ≤ l ≤ N . A derivada do segundo termo a direita, da equação (6.7), pode ser
aproximado por diferença finita, como segue.
u′ (0)x1−α h−α
≈ (u1 − u0 )
Γ (2 − α) Γ (2 − α)lα−1
Na aproximação do terceiro termo a direita, da equação (6.7), usamos a diferença finita central
de segunda ordem, em cada subintervalo, sendo

x x
u(2) (η) u(2) (x − η)
Z Z
1 1
α−1
dη = dη
Γ (2 − α) 0 (x − η) Γ (2 − α)0 (η)α−1
l−1 Z (j+1)h (2)
1 X u (x − η)
= dη
Γ (2 − α) j=0 jh (η)α−1
l−1
h−α X
(ul−j+1 − 2ul−j + ul−j−1 ) (j + 1)2−α − j 2−α

≈
Γ (3 − α) j=0

Desta forma, temos a aproximação para a derivada fracionária a direita (6.7) com 1 < α ≤ 2
dada por

h−α

GL α (1 − α)(2 − α)u0 (2 − α)(u1 − u0 )
D [u(x)][0,x] ≈ +
Γ (3 − α) lα lα−1
l−1
)
X
(ul−j+1 − 2ul−j + ul−j−1 ) (j + 1)2−α − j 2−α

+ (6.9)
j=0
A aproximação para a derivada fracionária a direita (6.8) com 1 < α ≤ 2 é obtida de forma
similar, sendo

80
6.2 Esquemas Numéricos

h−α

GL α (1 − α)(2 − α)uN (2 − α)(uN − uN −1 )
D [u(x)][x,L] ≈ α
+
Γ (3 − α) (N − l) (N − l)α−1
l−1
)
X
(ul+j−1 − 2ul+j + ul+j+1 ) (j + 1)2−α − j 2−α

+ (6.10)
j=0
Com isto temos que o termo de difusão pode ser discretizado substituindo (6.9) e (6.10) em (6.4).
De forma similar ao termo de difusão, vamos apresentar um esquema de aproximação com base
na Definição 6.1 e no Lema 6.1.1 para o termo de advecção. Para a aproximação do termo de advecção
vamos considerar as derivadas a esquerda e a direita de Grünwald-Letnikov, de ordem 0 < β < 1 no
intervalo [0, L].
Z x
GL α u(0)x−β 1 u′ (η)
D [u(x)][0,x] = + dη, (6.11)
Γ (1 − β) Γ (1 − β) 0 (x − η)β
Z L
GL α u(L)(L − x)−β 1 u′ (η)
D [u(x)][x,L] = + dη, (6.12)
Γ (1 − β) Γ (1 − β) x (η − x)β
Na aproximação do segundo termo a direita, da equação (6.11), usamos a diferença finita central
de segunda ordem, em cada subintervalo, sendo

x x
u′ (η) u′ (x − η)
Z Z
1 1
β
dη = dη
Γ (1 − β) 0 (x − η) Γ (1 − β)0 (η)β
l−1 Z (j+1)h ′
1 X u (x − η)
= dη
Γ (1 − β) j=0 jh (η)β
l−1
h−β X
(ul−j − ul−j−1 ) (j + 1)1−β − j 1−β

≈
Γ (2 − β) j=0

Obtemos assim, uma aproximação para a derivada fracionária a esquerda (6.11), com 0 < β < 1, é
dada por
( l−1
)
−β
GL α h (1 − β)u 0
X
(ul−j − ul−j−1 ) (j + 1)l−β − j l−β

D [u(x)][0,x] ≈ + (6.13)
Γ (2 − β) lβ j=0
Com procedimento similar, uma aproximação para a derivada fracionária a direita (6.12), com 1 <
β < 1, é dada por
( NX−l−1
)
−β
GL α h (1 − β)u N
(ul+j − ul+j+1 ) (j + 1)l−β − j l−β

D [u(x)][x,L] ≈ + (6.14)
Γ (2 − β) (N − l)β j=0
Com isto temos que o termo de difusão pode ser discretizado substituindo (6.13) e (6.14) em
(6.5).
Ao reunimos os resultados das aproximações dos termos de difusão e advecção, obtemos o

81
6.2 Esquemas Numéricos

sistema de equações diferenciais ordinárias na variável tempo dada por

Kα h−α

dul (1 − α)(2 − α)u0 (2 − α)(u1 − u0 )
≈ − +
2 cos πα

dt 2 Γ (3 − α)
lα lα−1
l−1
X (1 − α)(2 − α)uN
(ul−j+1 − 2ul−j + ul−j−1 ) (j + 1)2−α − j 2−α +

+
(N − l)α
j=0

−l−1
NX
(2 − α)(uN − uN −1 ) 2−α 2−α

+ + (u l+j−1 − 2u l+j + u l+j+1 ) (j + 1) − j
(N − l)α−1 
j=0

l−1
Kβ h−β  (1 − β)u
0
X h
1−β 1−β
i
− + (u l−j+1 − u l−j ) (j + 1) − j
2 cos πα

2 Γ (2 − α)
 lβ
j=0

−l−1
NX
(1 − β)uN h
1−β 1−β
i 
+ + (ul+j−1 − u l+j ) (j + 1) − j (6.15)
(N − l)β 
j=0

onde l = 1, 2, . . . , N − 1. Existem várias técnicas para solucionar o sistema de equações diferenciais

(6.15). Nos exemplo que iremos apresentar usamos a biblioteca Scify do python [5]. Aplicaremos o
método implı́cito BDF (Backward-Differentiation Formulas) da função solve_ivp .

6.2.2 Esquema de aproximação Grünwald deslocado

Proposto por Meerschaert e Tadjeran [29], o esquema de aproximação Grünwald deslocado para
a discretização da derivada fracionária bilateral corrige o problema de estabilidade apresentado pelo
esquema de Grünwald padrão na discretização do termo de difusão. Para o termo de difusão, vamos
considerar a discretização da derivada fracional de Riesz aplicando o esquema de Grünwald deslocado
l+1
GL α 1 X
D [u(x)][0,x] = α gj ul−j+1 + O(h) (6.16)
h j=0
N −l+1
GL α 1 X
D [u(x)][x,L] = α gj ul+j−1 + O(h) (6.17)
h j=0
(6.18)
onde os coeficientes gj são definidos por
α(α − 1) · · · (α − j + 1)
g0 = 1, gj = (−1)j
j!
com j = 1, 2, . . . , N .
Para o termo de advecção, a única mudança que temos que fazer é a variação dos ı́ndices em u.
Vamos considerar a discretização da derivada fracionária de Riesz aplicando o esquema de Grünwald

82
6.2 Esquemas Numéricos

padrão
l+1
GL α 1 X
D [u(x)][0,x] = α wj ul−j + O(h) (6.19)
h j=0
N −l+1
GL 1 X
Dα [u(x)][x,L] = wj ul+j + O(h) (6.20)
hα j=0
(6.21)
onde os coeficientes wj são definidos por
β(β − 1) · · · (β − j + 1)
w0 = 1, wj = (−1)j
j!
com j = 1, 2, . . . , N .
Com isto temos que o seguinte sistema de equações diferenciais ordinárias que será resolvido
pelo método implı́cito BDF (Backward Differentiation Formulas.) da função solve_ivp na variável
tempo
   
l+1 −l+1
NX l+1 −l+1
NX
dul Kα h−α X Kβ h−β X
≈− πα
 gj ul−j+1 + gj ul+j−1  −  wj ul−j + wj ul+j 
dt 2 cos 2 2 cos πβ
j=0 j=0 2 j=0 j=0
(6.22)

6.2.3 Esquema de aproximação por Transformação Matricial

Com base no trabalho de M. Ilić et al. [15], Yang e Tuner [52] propuseram um esquema
de transformação matricial para uma equação de difusão fracionária com condições de fronteiras
homogêneas. Vamos considerar a representação matricial de (6.1)
dU
= −ηα Aα/2 U − γβ Aβ/2 U (6.23)
dt
Kα Kβ L
onde ηα = α , γβ = β , h é o passo da discretização espacial definida por h = com
h h N

U = (u1 , u2 , . . . , uN −1 )T e A = tridiag(−1, 2, −1).

A matriz A é de ordem (N − 1) × (N − 1) e é simétrica definida positiva, ou seja, ela pode ser

decomposta como
A = PΛPT

onde Λ = diag(λ1 , λ2 , . . . , λN −1 ). λi , i = 1, . . . , N − 1 são os autovalores de A onde

kπ
λk = 2 − 2 cos , para k = 1, 2, . . . , n
n+1
As colunas da matriz P são os autovetores correspondentes e normalizados, onde
r
(j) 2 jkπ
vk = sin , para j = 1, 2, . . . , n
N +1 n+1

83
6.3 Exemplos Numéricos

Consequentemente a equação (6.23) se torna o seguinte sistema de equações diferenciais ordinárias

na variável tempo
dU
= −ηα PΛα/2 PT U − γβ PΛβ/2 PT U

(6.24)
dt
sujeita as condições iniciais
U(0) = (g1 , g2 , . . . , gN −1 )T

o qual pode ser resolvido pelo método implicito BDF.

6.3 Exemplos Numéricos

Nesta seção iremos apresentar exemplos numéricos para os esquemas apresentados na seção
anterior. Vamos considerar a equação de advecção-difusão fracionária to tipo de Riesz
Sejam 1 < α ≤ 2, 0 < β ≤ 1, 0 < t < 1 e 0 < x < π. A equação
∂u(x, t) ∂α ∂β
= Kα u(x, t) + K β u(x, t) (6.25)
∂t ∂|x|α ∂|x|β
sujeita as condições
u(0, t) = u(π, t) = 0 (6.26)
u(x, 0) = x2 (π − x), (6.27)
cuja a solução analı́tica dada por
+∞
8 4
− 3 sen(nπx)e−[Kα n +Kβ n ]t
α β
X
n+1
u(x, t) = 3
(−1) (6.28)
n=1
n n
Em todos os exemplos apresentados, o sistema de equações diferenciais ordinárias na variável
tempo será resolvido pelo método BDF, um método implı́cito de ordem variável (1 a 5) em múltiplas
etapas, baseado em uma fórmula de diferenciação reversa para a aproximação da derivada no tempo.
Na discretização dos termos de advecção e de difusão foi aplicado 50 pontos na variável espacial e
100 pontos na variável tempo. A ordem da derivada do termo de advecção usada é β = 0.4 e do termo
de difuão é α = 1.8. As constantes de advecção e difusão são Kα = Kβ = 0.25.

6.3.1 Aproximação L1-L2

Na Figura 6.1 temos a solução numérica obtida pelo esquema de aproximação L1-L2 apresentado
na Seção 6.2.1.
Fixando o tempo no nı́vel t = 0.51, o gráfico da Figura 6.2 apresenta a comparação da solução
numérica e solução analı́tica. No Apêndice D temos a implementação do esquema, usando a linguagem
Python. Aqui devemos observar que o erro apresentado na solução tem origem não só na discretização
dos termos de advecção e difusão, como também o erro na discretização da variável temporal.

84
6.3 Exemplos Numéricos

Figura 6.1: Aproximação L1L2

Figura 6.2: Comparação entre Aproximação L1-L2 com solução exata

6.3.2 Aproximação Grünwald Deslocado

A derivação fracionária de Grünwald apresenta problemas de aproximação quando 1 < α ≤ 2.
Isto faz com que o esquema de aproximação de Grünwald padrão (sem deslocamento) se torne
instável. Meerschaert e Tadjeran [29] propuseram o esquema de aproximação Grünwald deslocado
para a discretização da derivada fracionária bilateral, corrigindo o problema de estabilidade. Na
Figura 6.3 apresentamos a comparação entre a derivada de Grünwald padrão e a derivada de Grünwald
deslocada. As derivadas foram aplicadas em funções da forma gn (x) = sen(nπx), com n = 1, 2, 3, 4
e x ∈ [0, π]. Note que a solução analı́tica para a EADFR (6.28) é formada por esse tipo de função, o
que torna o esquema de aproximação padrão inadequado.
A solução numérica para o problema (6.25) sujeito as condições (6.27), obtida pelo esquema de
aproximação deslocado é apresentado na Figura 6.4. O código em Python encontra-se no Apêndice
E.
A Figura 6.5 apresenta a comparação entre a aproximação e a solução exata para o nı́vel de tempo
t = 0.51. .

6.3.3 Aproximação por Transformação Matricial

A estrutura do esquema de aproximação é tı́pica de matrizes tridiagonais simétricas com coefici-
entes constantes e aparece frequentemente em problemas de diferenças finitas (como a discretização
do operador Laplaciano unidimensional).

85
6.3 Exemplos Numéricos

(a) Derivada de sen(πx) (b) Derivada de sen(2πx)

(a) Derivada de sen(3πx) (b) Derivada de sen(4πx)

Figura 6.3: Comparação entre a derivada padrão e deslocada de Grünwald com α = 1.8

Figura 6.4: Aproximação Grünwald Deslocado

A solução numérica para o problema (6.25) sujeito as condições (6.27), obtida transformação
matricial é apresentado na Figura 6.6.
A comparação entre a aproximação e a solução exata para o nı́vel de tempo t = 0.51 é apresentada
na figura 6.7.
Neste capı́tulo apresentamos um tipo de equação diferencial parcial fracionária, onde as derivadas
fracionárias ocorriam na variável espacial. Para mais exemplos sugerimos o livro de B. Guo, X. Pu
e F. Huang [12]. Existem outras abordagem que devemos destacar. As derivadas fracionárias
são conhecidas por serem eficientes em fenômenos com memória, com as quais frequentemente
caracterizamos a dissolução de energia. A representação do modelo com operadores de derivadas
fracionárias, na variável tempo, demonstrou ser mais precisa e confiável quando comparada aos casos
de ordem inteira. Um exemplo é a equação de reação-difusão de Fisher não linear com derivada de
Caputo-Fabrizio fracionária no tempo [34].

86
6.3 Exemplos Numéricos

Figura 6.5: Comparação entre Aproximação Grünwald deslocado com solução exata

Figura 6.6: Aproximação por Transformação Matricial

Figura 6.7: Comparação entre Aproximação por Transformação Matricial com solução exata

87
Apêndice A Esquema Adams-Bashforth
O esquema de Adams-Bashforth é um esquema de passo múltiplo, onde os valores iniciais
y1 , . . . , yp−1 são obtidos por um esquema de passo simples. No código abaixo, temos a implementação,
onde os valores iniciais são obtidos pelo Esquema Numérico de Euler Melhorado e o número de passo
variando p = 1, 2, 3, 4. Isto garante ordem de precisão O(h3 ).

def AB1p(N, p, y0, t0, T,f):

"""
Método numérico AB1p (Adams-Bashforth de múltiplos passos)

Parâmetros:
N : número de pontos
p : ordem do segundo polinômio caracterı́stico
y0 : condição inicial
t0 : tempo inicial
T : tempo final
f : função f(t, y) para o cálculo numérico
EulerMelhorado : função que implementa o método de Euler Melhorado
"""

# Definição dos coeficientes beta

B = [Link]([
[2, 3, -1, 0, 0, 0, 0],
[12, 23, -16, 5, 0, 0, 0],
[24, 55, -59, 37, -9, 0, 0],
[720, 1901, -2774, 2616, -1274, 251, 0],
[1440, 4227, -7673, 9482, -698, 2627, -425]
])

# Definição de constantes
h = T / N

# Definição de vetores
y = [Link](N+1)
t = [Link](N+1)

t[0] = t0
y[0] = y0
# Método de Euler Melhorado
# Cálculo dos valores iniciais
for n in range(1, p+1):
t[n], y[n] = EulerMelhorado(t[n-1], y[n-1], h, f)

# Implementação do método AB1p

for m in range(p+1, N+1):
somaB = 0.0
for k in range(1, p+2):
somaB += B[p-1, k] * f(t[m-k], y[m-k])
y[m] = y[m-1] + h * somaB / B[p-1, 0]
t[m] = t[m-1] + h

return t, y

89
Apêndice B Esquema Adams-Moulton
Os Esquemas de Adams-Moulton, são esquemas implı́citos de passo p, ou seja o cálculo de ym
depende do cálculo de f (tm , ym ). Isto configura um esquema do tipo Preditor-Corretor, onde usamos
um esquema explı́cito de passo simples para calcular um primeiro valor de ym , que será “ corrigido”
pelo esquema de passo multiplo. Na implementação apresentada abaixo, usamos o Esquema de Euler
Melhorado para obter os valores iniciais e para atuar como esquema previsor.

def AM1p(N, p, y0, t0, T, f):

# Parâmetros:
# N: número de pontos
# p: ordem do segundo polinômio caracterı́stico
# y0: condição inicial
# t0: tempo inicial
# T: tempo final
# f: função externa f(x, y)

# Definição dos Coeficientes beta

B = [Link]([
[2, 1, 1, 0, 0, 0, 0],
[12, 5, 8, -1, 0, 0, 0],
[24, 9, 19, -5, 1, 0, 0],
[720, 251, 646, -264, 106, -19, 0],
[1440, 475, 1427, -798, 482, -173, 27]
])

# Definição de constantes
h = T / N

# Definição de vetores
y = [Link](N + 1)
t = [Link](N + 1)

t[0] = t0
y[0] = y0

# Método de Euler Melhorado para valores iniciais

for n in range(1, p+1):
t[n], y[n] = EulerMelhorado(t[n-1], y[n-1], h, f)
# Esquema de predição e correção
for m in range(p+1, N+1):
# Esquema Preditor
t[m], y[m] = EulerMelhorado(t[m-1], y[m-1], h, f)

# Esquema Corretor
somaB = 0.0
for k in range(1, p+2):
somaB += B[p, k] * f(t[m-k+1], y[m-k+1])

y[m] = y[m-1] + h * somaB / B[p, 0]

return t, y

91
Apêndice C Esquema de Adams-Bashforth-Moulton
Apresentamos uma versão do método de Adams-Bashforth-Moulton de um passo para resolver
um sistema de equação diferencial fracional para problemas na forma da integral Abel-Volterra. Na
integração utilizamos um esquema do tipo Newton-Côtes.

import numpy as np
from [Link] import gamma

# Definição da função f4t(x, y)

# O usuário deve definir essa função separadamente, pois não está no arquivo MATLAB forn
def f4t(yv, t):
# Transcrição direta da função f4t
# Retorna o valor negativo do primeiro elemento do vetor yv
f = -yv[0]

# Caso você queira experimentar as outras fórmulas, remova os comentários:

# f = -yv[0] - yv[3] + t**3 + 7 * t + 1 + 8 * t**(3/2) / [Link]([Link])
# f = t**(0.7) * ml(-t, 1.0, 1.7, 50) + [Link](-2 * t) - yv[0]**2 # Exemplo ex2

return f

def FDABM(N, alpha, y0, t0, T, S):

# Parâmetros:
# N: número de pontos
# alpha: ordem fracionária da EDOF
# y0: vetor com as condições iniciais
# t0: tempo inicial
# T: tempo final
# S: número de equações do sistema

# Definição de constantes
h = T / N
h_alpha = h ** alpha
m = int([Link](alpha))
gamma1 = gamma(alpha + 1)
gamma2 = gamma(alpha + 2)

# Definição de vetores
a = [Link](N) # Pesos Corretor
b = [Link](N) # Pesos Preditor
c = [Link](N) # Pesos Condições Iniciais Corretor
y = [Link]((N + 1, S))
t = [Link](N + 1)
fyt = [Link](N + 1)
ypred = [Link](S)
ym = [Link](S)

# Inicialização
y[0, :] = y0
t[0] = t0
fyt[0] = f4t(y[0, :], t[0])

# Cálculo dos pesos

for k in range(1, N):
a[k] = (k + 1) ** (alpha + 1) - 2 * k ** (alpha + 1) + (k - 1) ** (alpha + 1)
b[k] = (k ** alpha) - ((k - 1) ** alpha)
c[k] = (k - 1) ** (alpha + 1) - (k - 1 - alpha) * (k ** alpha)

# Loop principal
for j in range(1, N):
sumb = [Link](S)
suma = [Link](S)

for k in range(j):
for n in range(S - 1):
sumb[n] += b[j - k] * y[k + 1, n + 1]
sumb[S - 1] += b[j - k] * fyt[k + 1]

# Preditor
for k in range(S):
ypred[k] = y[0, k] + h_alpha * sumb[k] / gamma1

# Corretor
for k in range(j):
for n in range(S - 1):
suma[n] += a[j - k] * y[k + 1, n + 1]
suma[S - 1] += a[j - k] * fyt[k + 1]

93
t[j + 1] = (j + 1) * h

for k in range(S - 1):

ym[k] = c[j] * y[0, k + 1]
ym[S - 1] = c[j] * fyt[0]

for _ in range(3): # Aplica a correção três vezes

for n in range(S - 1):
ypred[n] = y[0, n] + h_alpha * (suma[n] + ym[n] + ypred[n + 1]) / gamma2
ypred[S - 1] = y[0, S - 1] + h_alpha * (suma[S - 1] + ym[S - 1] + f4t(ypred,

y[j + 1, :] = ypred
fyt[j + 1] = f4t(y[j + 1, :], t[j + 1])

return t, y

94
Apêndice D Esquema de aproximação L1-L2
Apresentamos uma versão do método de Aproximação L1-L2 descrito na Seção 6.2.1 aplicado na
Equação de Advecção-Difusão Fracionária do tipo de Riesz. Também apresentamos a implementação
da solução exata.

import numpy as np
from [Link] import gamma
from [Link] import solve_ivp
import [Link] as plt
from mpl_toolkits.mplot3d import Axes3D
from [Link] import quad
from [Link] import fft, ifft
import [Link] as plt

class RieszFractionalSolver:
def __init__(self, L, N, T, K_alpha, K_beta, alpha, beta):
"""
Inicializa o solver para equações fracionárias de Riesz

Parâmetros:
L : float - Comprimento do domı́nio espacial
N : int - Número de pontos da malha no espaço
T : float - Tempo final
K_alpha : float - Coeficiente de difusão fracionária
K_beta : float - Coeficiente de advecção fracionária
alpha : float - Ordem da derivada de difusao (1 < alpha <= 2)
beta : float - Ordem da derivada de advecção (0 < beta < 1)
"""
self.L = L
self.N = N
self.T = T
self.K_alpha = K_alpha
self.K_beta = K_beta
[Link] = alpha
[Link] = beta

# Discretização espacial
self.h = L / N
self.x = [Link](0, L, N+1)
# Coeficientes constantes
self.c_alpha = 1 / (2 * [Link]([Link] * alpha / 2))
self.c_beta = 1 / (2 * [Link]([Link] * beta / 2))

# Pré-calcular os pesos para L1 e L2

self._precompute_weights()

def _precompute_weights(self):
"""Pré-calcula os pesos para os esquemas L1 e L2"""
max_weights = self.N + 2

# Pesos para L2 (alpha em (1,2])

self.w_alpha = [Link](max_weights)
self.w_alpha[0] = 1
for j in range(1, max_weights):
self.w_alpha[j] = ((j+1)**([Link]) - j**([Link]))

# Pesos para L1 (beta em (0,1))

self.w_beta = [Link](max_weights)
self.w_beta[0] = 1
for j in range(1, max_weights):
self.w_beta[j] = ((j+1)**([Link]) - j**([Link]))

def _compute_L2_derivative(self, u):

"""Calcula a derivada fracionária de ordem alpha usando L2"""
N = self.N
h_alpha = self.h**[Link]
du = np.zeros_like(u)

for i in range(1, N): # Ignorar bordas (condições de Dirichlet)

# Termo à esquerda
left = 0
for j in range(i):
left += (u[i-j+1] - 2*u[i-j] + u[i-j-1]) * self.w_alpha[j]

# Termo à direita
right = 0
for j in range(N-i):
right += (u[i+j-1] - 2*u[i+j] + u[i+j+1]) * self.w_alpha[j]

96
# Combinação Riesz
du[i] = -self.c_alpha / gamma([Link]) * (left + right) / h_alpha

return du

def _compute_L1_derivative(self, u):

"""Calcula a derivada fracionária de ordem beta usando L1"""
N = self.N
h_beta = self.h**[Link]
du = np.zeros_like(u)

for i in range(1, N): # Ignorar bordas

# Termo à esquerda
left = 0
for j in range(i):
left += (u[i-j] - u[i-j-1]) * self.w_beta[j]

# Termo à direita
right = 0
for j in range(N-i):
right += (u[i+j] - u[i+j+1]) * self.w_beta[j]

# Combinação Riesz
du[i] = -self.c_beta / gamma([Link]) * (left + right) / h_beta

return du

def rhs(self, t, u):

"""Função do lado direito para o solver de ODEs"""
# Aplicar condições de contorno de Dirichlet
u_with_bc = [Link](self.N+1)
u_with_bc[1:self.N] = u
u_with_bc[0] = 0 # u(0,t) = 0
u_with_bc[self.N] = 0 # u(L,t) = 0

# Calcular os termos fracionários

d_alpha = self._compute_L2_derivative(u_with_bc)
d_beta = self._compute_L1_derivative(u_with_bc)

97
# Combinar os termos
return self.K_alpha * d_alpha[1:self.N] + self.K_beta * d_beta[1:self.N]

def solve(self, initial_condition):

"""
Resolve a equação fracionária

Parâmetros:
initial_condition : função - Função u(x,0) = g(x)

Retorna:
sol : Objeto solution de solve_ivp
"""
u0 = initial_condition(self.x[1:self.N])
sol = solve_ivp([Link], [0, self.T], u0, method=’BDF’,
t_eval=[Link](0, self.T, 100))
return sol

def plot_solution(x, t, u, alpha, beta):

"""Visualização 3D e 2D da solução"""
# Criar grid para 3D
T, X = [Link](t, x)
# Plot 3D
fig = [Link](figsize=(14, 6))
ax1 = fig.add_subplot(121, projection=’3d’)
surf = ax1.plot_surface(X, T, u, cmap=’viridis’, rstride=1, cstride=1,
linewidth=0, antialiased=True)
ax1.set_xlabel(’Posição x’, fontsize=12)
ax1.set_ylabel(’Tempo t’, fontsize=12)
ax1.set_zlabel(’u(x,t)’, fontsize=12)
ax1.set_title(f’Solução 3D alpha, beta)’, fontsize=14)
[Link](surf, ax=ax1, shrink=0.5, aspect=5)

# Plot 2D em tempos selecionados

ax2 = fig.add_subplot(122)
for i in range(0, len(t), max(1, len(t)//5)):
[Link](x,u[:, i], label=f’t = {t[i]:.2f}’)
ax2.set_xlabel(’x’, fontsize=12)
ax2.set_ylabel(’u(x,t)’, fontsize=12)

98
ax2.set_title(’ Evolução Temporal’, fontsize=14)
[Link]()
[Link](True)

plt.tight_layout()
[Link]()

class ExactFractionalSolution:
def __init__(self, L, alpha, beta, K_alpha, K_beta):
"""
Inicializa o solver para a solução exata

Parâmetros:
L : float - Comprimento do domı́nio
alpha : float - Ordem da derivada de difusão (1 < alpha <= 2)
beta : float - Ordem da derivada de advecção (0 < beta < 1)
K_alpha : float - Coeficiente de difusão
K_beta : float - Coeficiente de advecção
"""
self.L = L
[Link] = alpha
[Link] = beta
self.K_alpha = K_alpha
self.K_beta = K_beta

# Constantes
self.c_alpha = 1 / (2 * [Link]([Link] * alpha / 2)) if alpha != 1 else 1.0
self.c_beta = 1 / (2 * [Link]([Link] * beta / 2)) if beta != 1 else 1.0

def _eigenfunction(self, n, x):

"""Autofunção do operador de Riesz com condições de Dirichlet"""
return [Link](n * [Link] * x / self.L)

def _eigenvalue(self, n):

"""Autovalor para o modo n"""
return (n * [Link] / self.L) ** 2

def _compute_coefficient(self, n, initial_condition):

"""Calcula o coeficiente de Fourier b_n"""

99
integrand = lambda x: initial_condition(x) * self._eigenfunction(n, x)
b_n, _ = quad(integrand, 0, self.L)
return (2 / self.L) * b_n

def exact_solution(self, x, t, initial_condition, N_terms=50):

"""
Calcula a solução exata via expansão em série

Parâmetros:
x : array - Pontos espaciais onde calcular a solução
t : float - Tempo onde calcular a solução
initial_condition : função - u(x,0) = g(x)
N_terms : int - Número de termos na série

Retorna:
u : array - Solução em cada ponto x no tempo t
"""
u = np.zeros_like(x)

for n in range(1, N_terms + 1):

lambda_n = self._eigenvalue(n)
b_n = self._compute_coefficient(n, initial_condition)
decay = [Link](-self.K_alpha * lambda_n**([Link]/2) * t )
decay -= self.K_beta * lambda_n**([Link]/2) * t
u += b_n * self._eigenfunction(n, x) * decay

return u

def exact_solution_series(self, x, t_values, initial_condition, N_terms=50):

"""
Calcula a solução exata para múltiplos tempos

Parâmetros:
x : array - Pontos espaciais
t_values : array - Tempos para calcular
initial_condition : função - u(x,0) = g(x)
N_terms : int - Número de termos na série

Retorna:
u_matrix : array 2D - Solução em cada (x,t)

100
"""
u_matrix = [Link]((len(x),len(t_values)))

for i, t in enumerate(t_values):
u_matrix[:, i] = self.exact_solution(x, t, initial_condition, N_terms)

return u_matrix

101
Apêndice E Esquema de aproximação Grünwald
deslocado
Proposto por Meerschaert e Tadjeran [29], o esquema de aproximação Grünwald deslocado é
aplicado no termo difusivo, enquanto termo de advecção é aproximado pelo esquema padrão.

class ShiftedGrunwaldSolver:
def __init__(self, L, N, alpha, beta, K_alpha, K_beta):
"""
Inicializa o solver com deslocamento adaptativo para cada termo

Parâmetros:
L : float - Comprimento do domı́nio
N : int - Número de pontos da malha (excluindo bordas)
alpha : float - Ordem da derivada difusiva (1 < alpha <= 2)
beta : float - Ordem da derivada advectiva (0 < beta < 1)
K_alpha : float - Coeficiente de difusão
K_beta : float - Coeficiente de advecção
"""
self.L = L
self.N = N
[Link] = alpha
[Link] = beta
self.K_alpha = K_alpha
self.K_beta = K_beta
self.h = L / (N + 1) # Tamanho do passo espacial

# Pré-calcular coeficientes com deslocamentos distintos

self._compute_diffusion_coeffs() # Deslocamento 1 para termo difusivo
self._compute_advection_coeffs() # Sem deslocamento para termo advectivo

# Constantes de normalização
self.c_alpha = 1 / (2 * [Link]([Link] * alpha / 2)) if alpha != 1 else 1.0
self.c_beta = 1 / (2 * [Link]([Link] * beta / 2)) if beta != 1 else 1.0

def _compute_diffusion_coeffs(self):
"""Coeficientes para termo difusivo com deslocamento"""
max_k = self.N + 2
self.g_alpha = [Link](max_k)
self.g_alpha[0] = 1
for k in range(1, max_k):
self.g_alpha[k] = self.g_alpha[k-1] * (k - 1 - [Link] ) / k

def _compute_advection_coeffs(self):
"""Coeficientes para termo advectivo sem deslocamento"""
max_k = self.N + 2
self.g_beta = [Link](max_k)
self.g_beta[0] = 1
for k in range(1, max_k):
self.g_beta[k] = self.g_beta[k-1] * (k - 1 - [Link] ) / k

def _diffusive_term(self, u):

"""Termo difusivo com deslocamento delta=1"""
N = len(u) - 1
diff_term = np.zeros_like(u)

for i in range(1, N):

# Soma à esquerda (deslocada +1)
left = 0.0
for k in range(i + 1): # i+1 termos devido ao deslocamento
weight = self.g_alpha[k]
idx = i - k + 1 # Deslocamento positivo
left += weight * u[idx]

# Soma à direita (deslocada +1)

right = 0.0
for k in range(N - i + 1): # N-i+1 termos devido ao deslocamento
weight = self.g_alpha[k]
idx = i + k - 1 # Deslocamento negativo
right += weight * u[idx]

diff_term[i] = -self.c_alpha * (left + right) / (self.h ** [Link])

return diff_term

def _advective_term(self, u):

"""Termo advectivo sem deslocamento (delta=0)"""
N = len(u) - 1
adv_term = np.zeros_like(u)

103
for i in range(1, N):
# Soma à esquerda (não deslocada)
left = 0.0
for k in range(i): # i termos (sem deslocamento)
weight = self.g_beta[k]
idx = i - k # Sem deslocamento
left += weight * u[idx]

# Soma à direita (não deslocada)

right = 0.0
for k in range(N - i): # N-i termos (sem deslocamento)
weight = self.g_beta[k]
idx = i + k # Sem deslocamento
right += weight * u[idx]

adv_term[i] = -self.c_beta * (left + right) / (self.h ** [Link])

return adv_term

def rhs(self, t, u):

"""Função do lado direito para solve_ivp"""
# Aplicar condições de contorno (u[0] = u[-1] = 0)
u_full = [Link](self.N + 2)
u_full[1:-1] = u

# Calcular termos com deslocamentos distintos

diff_term = self.K_alpha * self._diffusive_term(u_full)
adv_term = self.K_beta * self._advective_term(u_full)

# Retornar apenas pontos internos

return diff_term[1:-1] + adv_term[1:-1]

def solve(self, initial_condition, t_span, t_eval):

"""Resolve a equação no intervalo temporal"""
# Malha espacial (pontos internos)
x = [Link](self.h, self.L - self.h, self.N)
u0 = initial_condition(x)

# Resolver sistema de EDOs

104
sol = solve_ivp([Link], t_span, u0, t_eval=t_eval, method=’BDF’)

# Retornar pontos espaciais (incluindo bordas) e solução

x_full = [Link](0, self.L, self.N + 2)
return x_full, sol

105
Apêndice F Esquema de aproximação por
Transformação Matricial
Yang e Tuner [52] propuseram um esquema de transformação matricial para uma equação de
difusão fracionária com condições de fronteiras homogêneas. Aplicamos este esquema para a equação
de advecção-difusão.

import numpy as np
from [Link] import solve_ivp
import [Link] as plt
from mpl_toolkits.mplot3d import Axes3D
from [Link] import gamma
from [Link] import quad
from [Link] import fft, ifft

class MTMFractionalSolver:
def __init__(self, L, N, alpha, beta, K_alpha, K_beta):
"""
Inicializa o solver MTM para equação de advecção-difusão fracionária

Parâmetros:
L : float - Comprimento do domı́nio
N : int - Número de pontos da malha (excluindo bordas)
alpha : float - Ordem da derivada de difusão (1 < \alpha <= 2)
beta : float - Ordem da derivada de advecção (0 < beta < 1)
K_alpha : float - Coeficiente de difusão
K_beta : float - Coeficiente de advecção
"""
self.L = L
self.N = N
[Link] = alpha
[Link] = beta
self.K_alpha = K_alpha
self.K_beta = K_beta
self.h = L / N

# Construir matrizes
self._build_matrices()
# Pré-computar as matrizes fracionárias
self._compute_fractional_matrices()

def _build_matrices(self):
"""Constrói as matrizes do Laplaciano e do gradiente discretos"""
# Matriz do Laplaciano (difusão)
main_diag = 2 * [Link](self.N)
off_diag = -1 * [Link](self.N - 1)
self.A = ([Link](main_diag) + [Link](off_diag, 1) + [Link](off_diag, -1))

def _compute_fractional_matrices(self):
"""Calcula as matrizes fracionárias via diagonalização"""
# Diagonalizar A = P D PˆT
self.eigvals_A = [Link](self.N)
for k in range(0,self.N):
self.eigvals_A[k] = 2 - 2*[Link]((k+1)*[Link]/(self.N+1))

[Link] = [Link]((self.N,self.N))
for j in range(0,self.N):
for k in range(0,self.N):
[Link][k][j] = [Link]((j+1)*(k+1)*[Link]/(self.N+1))

# Normalização
[Link] = [Link]/[Link]((self.N+1)/2)

# Calcular Aˆ{alpha/2} = P Dˆ{alpha/2} PˆT

Lambda_alpha = [Link](self.eigvals_A ** ([Link] / 2))
self.A_alpha = [Link] @ Lambda_alpha @ [Link].T

# Diagonalizar B para obter Bˆbeta

Lambda_beta = [Link](self.eigvals_A ** ([Link] / 2))
self.B_beta = [Link] @ Lambda_beta @ [Link].T

def rhs(self, t, u):

"""Função do lado direito para solve_ivp"""
return -self.K_alpha /(self.h ** alpha) * (self.A_alpha @ u) - self.K_beta /(sel

107
def solve(self, initial_condition, t_span, t_eval):
"""
Resolve a equação fracionária

Parâmetros:
initial_condition : função - u(x,0) = g(x)
t_span : tuple - (t0, tf)
t_eval : array - Tempos onde a solução é calculada

Retorna:
x : array - Pontos da malha espacial (incluindo bordas)
sol : Objeto solution de solve_ivp
"""
# Discretizar condição inicial (pontos internos)
x_internal = [Link](self.h, self.L - self.h, self.N)
u0 = initial_condition(x_internal)

# Resolver o sistema de EDOs

sol = solve_ivp([Link], t_span, u0, t_eval=t_eval, method=’LSODA’)

# Adicionar condições de contorno (zeros nas bordas)

x_full = [Link](0, self.L, self.N + 2)
return x_full, sol

def plot_comparison(x, t, u_mtm, alpha, beta):

"""Comparação visual entre MTM e Grünwald"""
fig = [Link](figsize=(14, 6))

# Plot 3D para MTM

ax1 = fig.add_subplot(121, projection=’3d’)
T, X = [Link](t, x)
ax1.plot_surface(X, T, u_mtm, cmap=’viridis’, alpha=0.8)
ax1.set_title(f’MTM (alpha={alpha}, beta ={beta})’)
ax1.set_xlabel(’x’)
plt.tight_layout()
[Link]()

def plot_solution(x, t, u, alpha, beta):

"""Visualização 3D e 2D da solução"""
# Criar grid para 3D

108
T, X = [Link](t, x)

# Plot 3D
fig = [Link](figsize=(14, 6))
ax1 = fig.add_subplot(121, projection=’3d’)
surf = ax1.plot_surface(X, T, u, cmap=’viridis’, rstride=1, cstride=1,
linewidth=0, antialiased=True)
ax1.set_xlabel(’Posição x’, fontsize=12)
ax1.set_ylabel(’Tempo t’, fontsize=12)
ax1.set_zlabel(’u(x,t)’, fontsize=12)
ax1.set_title(f’Solução 3D (alpha={alpha}, beta={beta})’, fontsize=14)
[Link](surf, ax=ax1, shrink=0.5, aspect=5)

# Plot 2D em tempos selecionados

ax2 = fig.add_subplot(122)
for i in range(0, len(t), max(1, len(t)//5)):
[Link](x, u[:, i], label=f’t = {t[i]:.2f}’)
ax2.set_xlabel(’x’, fontsize=12)
ax2.set_ylabel(’u(x,t)’, fontsize=12)
ax2.set_title(’Evolução Temporal’, fontsize=14)
[Link]()
[Link](True)

plt.tight_layout()
[Link]()

109
Bibliografia
[1] D. Baleanu et al. Fractional Calculus: Models and Numerical Methods. Series on Complexity,
Nonlinearity and Chaos. Singapore: World Scientific, 2017.
[2] M. Caputo. ≪Linear models of dissipation whose Q is almost frequency independent II≫. Em:
Geophys J. Royal Astr. Soc. 13 (1967), pp. 529–539.
[3] B. C. Carlson. Special functions of applied mathematics. New York: Academic Press, 1977.
[4] A. Cayley. ≪Note on Riemann’s paper ”Versuch einer allgemeinen Auflösung der Integration
un Differentiation”≫. Em: Math. Ann. 16 (1880), pp. 81–82.
[5] SciPy community. SciPy Reference Guide. Release 1.8.1. Disponı́vel em https : / / docs .
[Link]/doc/scipy-1.8.1/[Link]. [Link]. Mai. de 2022.
[6] ≪The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using
Differential Operators of Caputo Type≫. Em: Lecture notes in mathematics. Ed. por K. Diethelm.
Springer Verlag, 2010.
[7] A. Erdélyi. ≪On fractional integration and its application on the theory of Hankel tran-
forms.≫ Em: Quart. J. Math. 11.44 (1940), pp. 293–303.
[8] A. Erdélyi e H. Kober. ≪Some remarks on Hankel tranforms≫. Em: Quart. J. Math. 11.43
(1940), pp. 212–221.
[9] L. Euler. ≪De progressionibus transcendibus seu quarum termini generales algebraicae dari
nequent≫. Em: Comm. Acad. Sci. Petropolitanae 5 (1738). Transleted to englis by S.G. Langton,
University of San Diego, [Link]/ langton, pp. 36–57.
[10] A. Kilbas G. Samko e O. Marichev. Fractional Integrals and derivatives: Teory and Applicati-
ons. Yverdon: Gordon e Breach Science Publishers, 1993.
[11] A. K. Grünwald. ≪Ueber ”begrenzte” Derivationen und deren Anwendung Z. angew≫. Em:
[Link] Phys. 12 (1867), pp. 441–480.
[12] B. Guo, X. Pu e F. Huang. Fracional Partial Differential Equations and Their Numerical
Solutions. Singapore: Word Scientific Publishing Co., 2015.
[13] O. Heaviside. Electrical papers. Vol. 8. London: The Macmillan Company, 1892.
[14] P. Henrici. Discrete variable methods in ordinary differential equations. Ed. por Camb. Univ.
Pr. New York: Encyclopedia of mathematics e its application., 1968.
[15] M. Ilić et al. ≪Numerical approximation of a fractional-in-space diffusion equation≫. Em: Fract.
Calculus Appl. Anal. 8.3 (2005), pp. 323–341.
[16] J. L. Lagrange. Sur une: Nouvelle espèce de calcul relatif. A la différentiation et a lı́ntégration
des quantites variable. Paris: Gauthier-Villars, 1849.
BIBLIOGRAFIA

[17] P. S. Laplace. Théorie analytique des probabilités. Paris: Imprimeur-Libraire pour les Mathématiques,
1820.
[18] H. Laurent. ≪Sur le calcul des dérivées à indicies quelconques.≫ Em: Nouv. Annales de
Mathématiques 3 (1884), pp. 240–252.
[19] A. V. Letnikov. ≪An explanation of the concepts of the theory of differentiation of arbitrary
index≫. Em: Moscow Matem. Sbornik 6 (1872), pp. 413–445.
[20] A. V. Letnikov. ≪Theory of differentiation with an arbitrary index≫. Em: Matem. Sbornik 3
(1868), pp. 1–66.
[21] J. Liouville. ≪Mémoire sur le théoréme des fonctions complémentaires.≫ Em: J. für reine und
angew. Math. 11 (1834), pp. 01–19.
[22] J. Liouville. ≪Mémoire sur quelques questions de géométrie et de mécanique, et sur um nouveau
genre de calcul pour résoudre ces questions≫. Em: J. l’Ecole Roy. Polytéchn 13 (1832), pp. 01–
69.
[23] E. R. Love. ≪Changing the order of integration≫. Em: J. Austral. Math. Soc. 9 (1970), pp. 421–
432.
[24] C. Lubich. ≪A stability analysis of convolution quadratures for Abel-Volterra integral≫. Em:
IMA J. Numer. Anal 6 (1) (1985), pp. 87–101.
[25] C. Lubich. ≪Discretized fractional calculus≫. Em: SIAM J. Math. Anal. 17.3 (1986), pp. 704–
719.
[26] C. Lubich. ≪Fractional linear multistep methods for Abel-Volterra integral equations of the
second kind.≫ Em: Math. Comp 45 (1985), pp. 463–469.
[27] C. Lubich. ≪On the stability of linear multistep methods for Volterra convolution equations≫.
Em: IMA J. Num. Anal. 3.4 (1983), pp. 439–465.
[28] A. Marchaud. ≪Sur les dérivées et sur les différences des fonctions des variables réelles.≫ Em:
J. Math. Pures et Appl. 6.4 (1927), pp. 371–382.
[29] M. M. Meerschaert e C. Tadjeran. ≪Finite difference approximations for two-sided spaced-
fractional partial differential equations≫. Em: Appl. Numer. Math. 56.1 (2006), pp. 80–90.
[30] K. S. Miller e B. Ross. An introduction to de fractionalcalculus and fractional diferential
equation. New York: Wiley-Interscience Publication, 1993.
[31] R.K. Miller. Nonlinear Volterra Integral Equations. Mathematics lecture note series. W. A.
Benjamin, 1971. isbn: 9780805370805.
[32] Feliz Manuel Barrão Minhós. Equações Diferenciais Ordinárias: Relatório sobre Unidade
Curricular. Rel. téc. Tech. Report. 2009.
[33] L. O’Shaughnessy. ≪Probem # 433≫. Em: Amer. Math. Monthly 25 (1918), pp. 172–173.
[34] K. M. Owolabi e A. Atangana. Numerical Methods for Fractional Differentiation. Springer
Series in Computational Mathematics. Singapore: Springer, 2019.

111
BIBLIOGRAFIA

[35] I. Podlubny. ≪Matrix approach to discrete fractional calculus.≫ Em: Fract. Calc. Appl. Anal 4
(2000), pp. 359–386.
[36] E. L. Post. ≪Discussion of the solution of (d/dx)1/2 y = y/x (Probem # 433)≫. Em: Amer.
Math. Monthly 26 (1919), pp. 37–39.
[37] R. L. Pouso. A simple proof of the Fundamental Theorem of Calculus for the Lebesgue integral.
Rel. téc. Tech. Report. 2012.
[38] M. Riesz. ≪L’ intégrales de Riemann-Liouville et le probléme de Cauchy.≫ Em: Acta Math. 81
(1936), pp. 1–223.
[39] M. Riesz. ≪L’ intégrales de Riemann-Liouville et potentiels≫. Em: Acta Litt. Acad. Sci. 9 (1938),
pp. 1–42.
[40] M. Riesz. ≪L’ intégrales de Riemann-Liouville et solution invariante du probléme de Cauchy
pour l’equation des ondes≫. Em: C. R. Congrés Intern. Math. 2 (1936), pp. 44–45.
[41] M. Riesz. ≪Potentiels de divers ordres et leurs fonctions de Green≫. Em: C. R. Congrés Intern.
Math. 2 (1936), pp. 62–63.
[42] ≪Fractional calculus and its application: proceedings of international conference held at Uni-
versity of New Haven≫. Em: Lecture notes in mathematics. Ed. por B. Ross. Springer Verlag,
1975.
[43] B. Ross. ≪Origins of fractional calculus and some applications≫. Em: Internat. J. Math. Statist.
Sci. 1.1 (1992), pp. 21–34.
[44] B. Ross. ≪The development of fractional calculus 1695-1900≫. Em: Historia Math 4 (1977),
pp. 75–89.
[45] O. Marichev S. Samko A. Kilbas. Fractional Integrals and derivatives: Theory and Applications.
Yverdon: Gordon e Breach Science Publishers, 1993.
[46] N. Y. Sonin. ≪On differentiation with arbitrary index.≫ Em: Moscow Matem. Sbornik 6 (1869),
pp. 1–38.
[47] K. Spindler. ≪A short proof of the formula of Faà di Bruno≫. Em: Elem. Math. 60 (2005),
pp. 33–35.
[48] G. Strang. Calculus. https : / / ocw . mit . edu / ans7870 / resources / Strang / Edited /
Calculus/[Link]. Massachusetts: Wellesley College, 1991.
[49] P.J. Torvik e R.L. Bagley. ≪On the appearance of the fractional derivativa in the behavior of
real materials≫. Em: [Link] Mechanics 51 (2) (1984), pp. 294–298.
[50] Y. Watanabe. ≪Notes on the generalized derivative of Riemann-Liouville and its applications
to Leibnitz’s formula, I and II≫. Em: Tóhoku Math. J. 34 (1931), pp. 8–41.
[51] H. Weyl. ≪Bemerkungen zum Begriff des Differentialquotienten gebrochener Ordnung.≫ Em:
Vierteljshr. Naturforsch. Gesellsch. Zürich 62 (1917), pp. 296–302.

112
BIBLIOGRAFIA

[52] Q. Yang, F. Liu e I. Turner. ≪Numerical methods for fractional partial differential equations
with Riesz space fractional derivatives≫. Em: Applied Mathematical Modelling 34 (2010),
pp. 200–218.

113
O Cálculo Fracionário surgue como uma generalização do cálculo classico de ordem inteira,
proposto por Newton e Leibniz, no século XVII. O tema começou a despertar interesse da comunidade
cientı́fica, em 1974, com a realização da primeira conferência internacional sobre Cálculo Fracionário,
promovida por Bertram Ross em New Havem, USA. Na década de 90 ocorre uma explosão de
publicações de livros e artigos sobre o tema, consolidando o Cálculo Fracionário como uma importante
ferramenta na análise de fenômenos, demonstrando ter propriedades e caracterı́sticas superiores ao
cálculo integro-deferencial clássico.
Neste trabalho fazemos uma breve revisão histórica da teoria do Cálculo Fracionário e suas
aplicações. São apresentados alguns resultados do Cálculo Fracionário, apresentando alguns exemplos
de equações diferenciais ordinárias fracionárias no contexto analı́tico e numérico.