Resolvendo numericamente equações diferenciais fracionárias

Solving fractional differential equations numerically

Resumos

Apresentamos neste trabalho uma breve revisão sobre o cálculo fracionário e também o método de Adam-Bashforth fracionário para solução de equações diferenciais fracionárias. O método de Adam-Bashforth fracionário é devido a Diethelm et al. [¹] e caracteriza-se pela sua eficiência, já que considera os efeitos de memória das derivadas fracionárias. A apresentação do tema é realizada de forma pedagógica, proporcionando a compreensão daqueles que são incipientes no assunto. Nesse contexto, exemplos são apresentados e algumas aplicações são discutidas.

Palavras-chaves:
Cálculo Fracionário; Derivada de Caputo; Método de Adam-Bashforth Fracionário

In this paper we present a brief review about fractional calculus and also the fractional Adam-Bashforth method for solving fractional differential equations. The fractional Adam-Bashforth method is due to Diethelm et al. [¹] and is characterized by its efficiency, as it considers the memory effects of fractional differentiation. The presentation is performed in a pedagogical perspective, delivered to beginner in the subject. In this context, some examples are presented and applications are discussed.

Keywords
Fractional Calculus; Caputo’s Derivative; Fractional Adam-Bashforth Method

1. Introdução

Embora pouco popular, o Cálculo Fracionário é mais antigo do que imaginamos. A sua origem é contemporânea do cálculo diferencial de ordem inteira, tradicionalmente estudado nos cursos universitários. O indício do surgimento dessa generalização do cálculo de ordem inteira é datado no século dezessete numa carta que L’Hôspital enviou a Leibniz indagando sobre a possibilidade de uma derivada de ordem $1/2$ [²]. A resposta dada por Leibniz culminou na introdução das primeiras definições de derivada e integral de ordem não-inteira. Em 1812, Laplace estabeleceu a primeira definição de derivada fracionária. A partir desse fato, eminentes matemáticos passaram a investigar as consequências da definição de derivadas fracionárias e aplicá-las na abordagem de problemas físicos. Nesse caminho, a primeira vez que a denominação “derivada fracionária” apareceu na literatura foi em 1819 nos trabalhos de Lacroix. No bojo das aplicações, um exemplo foi o matemático Abel, que em 1823 usou a derivada fracionária para estudar o conhecido problema da tautócrona. O estudo de Abel é considerado como a primeira aplicação prática do cálculo fracionário. Já em 1892, surge a definição de derivada e integral fracionária proposta por Riemann-Liouville. Tal definição é muito estudada até hoje, e tem sido aplicada na modelagem de diversos problemas. Contudo, a derivada fracionária de Riemann-Liouville mostrou-se um pouco obsoleta por alguns motivos, sendo um deles o fato da possibilidade da derivada fracionária de uma constante ser diferente de zero [³]. A partir desse problema, em 1967 emerge a derivada fracionária proposta por Michele Caputo [⁴], a qual soluciona parte dos desconfortos oriundos da versão de Riemann-Liouville, e principalmente por esse motivo, tem sido aplicada de forma abundante em problemas da física e da engenharia [³, ⁵, ⁶]. A ideia inicial de Caputo era desenvolver uma abordagem da viscoelasticidade usando cálculo fracionário, tendo muito sucesso em sua proposta. Uma característica do cálculo fracionário é que a derivada fracionária pode ser definida de muitas formas distintas, e geralmente, o resultado de um definição geralmente não coincide com o de outras. As definições apresentadas por Riemann-Liouville e Caputo são as mais conhecidas e serão abordadas neste trabalho, com ênfase dada na formulação de Caputo devido a proposta deste trabalho. O cálculo fracionário tem se desenvolvido bastante a partir de 1967 e tem sido utilizado para modelar diversas situações nas mais variadas áreas, quais sejam: física [⁷]; ciências mecânicas [⁸, ⁹, ¹⁰]; modelos matemáticos para ciências biológicas [¹¹, ¹²]; dentre muitas outras. As principais características da derivada fracionária que ampliam o seu horizonte de aplicações são a não-localidade e o efeito de memória [³]. Nesse sentido, um dos campos de estudo teórico do cálculo fracionário são as equações diferenciais fracionárias. Relevantes avanços já foram alcançados com a demonstração de teoremas e o desenvolvimento de técnicas de solução [⁵]; porém, há muito ainda por desenvolver sobretudo por haver equações diferenciais fracionárias relativamente simples para as quais não são conhecidas soluções analíticas. É justamente nesse ponto que surge a relevância deste artigo. Apresentaremos de forma pedagógica alguns elementos do cálculo fracionário e como principal constructo, abordaremos alguns métodos numéricos de solução de equações diferenciais fracionárias. Assim sendo, na seção 2 apresentamos os elementos matemáticos fundamentais para a compreensão do trabalho; além da definição de derivada fracionária segundo Caputo. Na seção 3 introduziremos as equações diferenciais fracionárias. Na seção 4 apresentaremos um método numérico para equações fracionárias e análogo ao método de Adams-Bashforth para equações diferenciais de ordem inteira. Três aplicações de áreas distintas serão apresentadas na seção 5. Enquanto na seção 6 elencamos as considerações finais e perspectivas.

2. Preliminares

Para que a derivada fracionária possa ser definida e aplicada será necessário definir e explorar algumas funções especiais que ajudarão a construir o conhecimento necessário para o desenvolvimento dessa. A seguir revisaremos a Funções Gama, a Função Beta e a Função de Mittag-Leffler.

2.1. Função Gama

A Função Gama, ou integral de Euler de segunda espécie, é uma das funções especiais mais importantes na matemática, tendo sido apresentada pela primeira vez pelo matemático Euler em 1730, como resultado de uma pesquisa sobre uma forma de interpolação do fatorial de um número [²]. Foi estudada por diversos outros matemáticos, mas, Adrian Marie Legendre, no ano de 1809, quem a definiu na forma que trabalhamos nos dias atuais. Para $z\in ℂ$, com $\mathrm{\Re }(z)>0$, a Função Gama é definida por [¹³]

(1) $$\mathrm{\Gamma }(z)={\int }_0^{\mathrm{\infty }}{e}^{-t}{t}^{z-1}\mathit{d}t.$$

Um dos fatos notáveis da Função Gama é que ela generaliza o fatorial de um número natural n. É imediato verificar a partir de uma integração por partes que

(2) $$\mathrm{\Gamma }(z+1)=z\mathrm{\Gamma }(z)$$

para qualquer z. Consequentemente, por indução, para $n\in \mathbb{N}$ temos que

(3) $$\mathrm{\Gamma }(n+1)=n!,n\in \mathbb{N}.$$

2.2. Função Beta

A Função Beta foi introduzida pela primeira vez por Euler e também é conhecida por integral de Euler de primeira espécie, definida como [¹³]

(4) $$\beta (q,p)={\int }_0^1{t}^{p-1}{(1-t)}^{q-1}\mathit{d}t,$$

em que $\mathrm{\Re }(q),\mathrm{\Re }(p)>0$. A Função Beta apresenta diversas propriedades interessantes as quais podemos destacar: (i) a comutatividade

$$\beta (q,p)=\beta (p,q),$$

cuja demonstração é imediata, e (ii) a relação entre a Função Beta e a Função Gama

(5) $$\beta (q,p)=\frac{\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)}{\mathrm{\Gamma }(q+p)}.$$

Para demonstrarmos essa identidade consideramos o produto

(6) $$\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)={\int }_0^{\mathrm{\infty }}{e}^{-t}{t}^{q-1}\mathit{d}t{\int }_0^{\mathrm{\infty }}{e}^{-z}{z}^{p-1}\mathit{d}z$$

Introduzindo as mudanças de variáveis $t=u^2$ e $z=v^2$, na Eq. (6), pode-se escrever,

$$\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)=4{\int }_0^{\mathrm{\infty }}{e}^{-u^2}{u}^{2q-1}\mathit{d}u{\int }_0^{\mathrm{\infty }}{e}^{-v^2}{v}^{2p-1}\mathit{d}v$$ $$=4{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{e}^{-(u^2+v^2)}{u}^{2q-1}{v}^{2p-1}\mathit{d}u\mathit{d}v.$$

Usando as coordenadas polares no plano $u=r\cos \theta$ e $v=r\mathrm{sen}\theta$, em que o Jacobiano da transformação é r, então $dudv=rdrd\theta$ e

$$\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)=\left(2{\int }_0^{\mathrm{\infty }}{e}^{-r^2}{r}^{2p+2q-1}\mathit{d}r\right)$$ $$\times \left(2{\int }_0^{\pi /2}{(\cos \theta )}^{2p-1}{(\mathrm{sen}\theta )}^{2q-1}\mathit{d}\theta \right).$$

Aqui vamos notar que

$$\beta (q,p)=2{\int }_0^{\pi /2}{(\mathrm{sen}\theta )}^{2p-1}{(\cos \theta )}^{2q-1}\mathit{d}\theta .$$

Pois, fazendo a mudança de $t={\cos }^2\theta$ na Eq. (4), segue que,

$$\beta (q,p)=-{\int }_{\pi /2}^0[{(\cos \theta )}^{2p-2}{(\mathrm{sen}\theta )}^{2q-2}\times$$ $$\times 2\mathrm{s}\mathrm{e}\mathrm{n}\theta \cos \theta d\theta ],$$

em que ${(\mathrm{sen}\theta )}^{2q-2}=1-t$ e $2\mathrm{s}\mathrm{e}\mathrm{n}\theta \cos \theta d\theta =dt$. Assim, temos que

$$\beta (q,p)=2{\int }_0^{\pi /2}{(sen\theta )}^{2p-1}{(cos\theta )}^{2q-1}\mathit{d}\theta .$$

Portanto,

(7) $$\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)=2\beta (q,p){\int }_0^{\mathrm{\infty }}{e}^{-r^2}{r}^{2p+2q-1}\mathit{d}r.$$

Finalmente, fazendo a mudança de variáveis $r^2=y$ e 2rdr = dy em Eq. (7) temos

$$\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)=\beta (q,p){\int }_0^{\mathrm{\infty }}{e}^{-y}{y}^{p+q-1}\mathit{d}y.$$

Logo, chegamos à identidade desejada

$$\beta (q,p)=\frac{\mathrm{\Gamma }(q)\mathrm{\Gamma }(p)}{\mathrm{\Gamma }(q+p)}.$$

2.3. Função de Mittag-Leffler

A Função de Mittag-Leffler é conhecida como uma generalização da função exponencial. Essa função, conforme introduzida por Leffler, é uma função de variável complexa que depende de um parâmetro complexo $\alpha$, onde $\mathrm{\Re }(\alpha )>0$, e é definida na forma de uma série de potências [¹⁴]

(8) $$E_{\alpha }(t)=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{\mathrm{\Gamma }(\alpha k+1)}.$$

Note que no caso em que $\alpha =1$, tem-se:

$$E_1(t)=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{\mathrm{\Gamma }(k+1)}=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{t!}=e^t,$$

o que permite-nos observar que a Função de Mittag-Leffler admite, como caso particular, a função exponencial.

2.3.1. Função de Mittag-Leffler de dois parâmetros

A Função de Mittag-Leffler de dois parâmetros, $E_{\alpha ,\beta }(t)$, foi introduzida por Wiman [¹⁴] e é uma função de variável complexa que depende de dois parâmetros complexos, $\alpha$ e $\beta$, onde $\mathrm{\Re }(\alpha ),\mathrm{\Re }(\beta )>0$, dada pela equação a seguir,

(9) $$E_{\alpha ,\beta }(t)=\sum _{k=0}^{\mathrm{\infty }}\frac{t^k}{\mathrm{\Gamma }(\alpha k+\beta )}.$$

A condição $\mathrm{\Re }(\alpha ),\mathrm{\Re }(\beta )>0$ é exigida para garantir a convergência da série dada na Eq. (9).

É imediato observar que quando $\beta =1$ a Eq. (9) se reduz a Função de Mittag-Leffler de um parâmetro, ou seja,

$$E_{\alpha \mathrm{,1}}(t)=E_{\alpha }(t).$$

2.3.2. Transformada de Laplace da Função deMittag-Leffler

Nas seções posteriores necessitaremos da transformada de Laplace da Função de Mittag-Leffler. Por essa razão, apresentamos neste momento os passos de sua obtenção. A transformada de Laplace de uma função $f(t)$, definida para $t\ge 0$, é dada pela integral

$$ℒ\{f(t)\}={\int }_0^{\mathrm{\infty }}{e}^{-st}f(t)\mathit{d}t,$$

desde que a integral convirja. Com isso, podemos calcular a transformada de Laplace da seguinte generalização da Função de Mittag-Leffler com dois parâmetros

$$ℒ\left\{t^{\beta -1}{E}_{\alpha ,\beta }(\lambda t^{\alpha })\right\}={\int }_0^{\mathrm{\infty }}{e}^{-st}\sum _{k=0}^{\mathrm{\infty }}\frac{{\lambda }^k{t}^{\alpha k+\beta -1}}{\mathrm{\Gamma }(\alpha k+\beta )}ds,$$

em que $\lambda$ é um parâmetro qualquer.

Como a Transformada de Laplace de $t^n$ é $ℒ\{t^n\}=\frac{n!}{s^{n+1}}$, obtemos

$$ℒ\left\{t^{\beta -1}{E}_{\alpha ,\beta }(\lambda t^{\alpha })\right\}=\sum _{k=0}^{\mathrm{\infty }}\frac{{\lambda }^k}{\mathrm{\Gamma }(\alpha k+\beta )}{\int }_0^{\mathrm{\infty }}{e}^{-st}{t}^{\alpha k+\beta -1}\mathit{d}t,$$ $$=\sum _{k=0}^{\mathrm{\infty }}\frac{{\lambda }^k}{\mathrm{\Gamma }(\alpha k+\beta )}\frac{\mathrm{\Gamma }(\alpha k+\beta )}{s^{\alpha k+\beta }}.$$

Logo,

$$ℒ\left\{t^{\beta -1}{E}_{\alpha ,\beta }(\lambda t^{\alpha })\right\}=\sum _{k=0}^{\mathrm{\infty }}\frac{{\lambda }^k}{s^{\alpha k+\beta }}$$ $$=\frac{1}{s^{\beta }}\sum _{k=0}^{\mathrm{\infty }}{\left(\frac{\lambda }{s^{\alpha }}\right)}^k.$$

Se recordarmos da série geométrica que é dada por

$$\frac{1}{1-x}=\sum _{k=0}^{\mathrm{\infty }}{x}^k,$$

quando $|x|<1$, temos que

(10) $$ℒ\left\{t^{\beta -1}{E}_{\alpha ,\beta }(\lambda t^{\alpha })\right\}=\frac{1}{s^{\beta }}\frac{1}{1-\lambda /s^{\alpha }},$$

contanto que $|\lambda /s^{\alpha }|<1.$

Finalmente, reescrevendo a Eq. (10), temos que a transformada de Laplace da Função de Mittag-Leffler é dada por,

(11) $$ℒ\left\{t^{\beta -1}{E}_{\alpha ,\beta }(\lambda t^{\alpha })\right\}=\frac{s^{\alpha -\beta }}{s^{\alpha }-\lambda }.$$

Este resultado será relevante na sequência do trabalho.

2.4. A derivada fracionária de Caputo

Nesta seção apresentaremos a definição de derivada fracionária segundo Caputo. Dedicaremos a nossa atenção a este tipo de derivada fracionária, em detrimento das outras possibilidades, devido a sua importância na análise de sistemas físicos [³].

Definição: A derivada fracionária de Caputo de ordem $\alpha$ é definida por [³]

(12) $${}_b{}^C\mathbf{D}_t^{\alpha }f(t)=\frac{1}{\mathrm{\Gamma }(n-\alpha )}{\int }_b^t\frac{f^{(n)}(u)}{{(t-u)}^{\alpha -n+1}}\mathit{d}u,$$

em que $\alpha \in {ℝ}_{>0},$ $b\in ℝ$, $n=\lceil \alpha \rceil$ e $f^{(n)}(u)=d^{n}f(u)/d{u}^n$ é a n-ésima derivada de ordem inteira. As derivadas de Caputo são válidas para $\alpha \in ℂ$, desde que $\mathrm{\Re }(\alpha )>0$, no entanto, aqui será considerado apenas $\alpha \in {ℝ}_{>0}$. Daqui em diante, por não haver confusão com outra definição de derivada fracionária, escreveremos o operador derivada de Caputo na forma ${}_0^C{\mathbf{D}}_t^{\alpha }={\mathbf{D}}_t^{\alpha }.$

Para fins didáticos, vamos calcular a derivada fracionária segundo Caputo da função potência $t^r$, com $\alpha \in {ℝ}_{>0}$. Por definição ela é dada por

$${\mathbf{D}}_t^{\alpha }[t^r]=\frac{1}{\mathrm{\Gamma }(n-\alpha )}{\int }_0^t\frac{{[u^r]}^{(n)}}{{(t-u)}^{\alpha -n+1}}\mathit{d}u,$$

em que $n=\lceil \alpha \rceil .$ Então,

$${\mathbf{D}}_t^{\alpha }[t^r]=\frac{1}{\mathrm{\Gamma }(n-\alpha )}{\int }_0^t\frac{r(r-1)\mathrm{\dots }(r-n+1)u^{r-n}}{{(t-u)}^{\alpha -n+1}}\mathit{d}u$$ $$=\frac{\mathrm{\Gamma }(r+1)}{\mathrm{\Gamma }(n-\alpha )\mathrm{\Gamma }(r-n+1)}$$ $$\times {\int }_0^t{u}^{r-n}{(t-u)}^{n-\alpha -1}du$$ $$=\frac{\mathrm{\Gamma }(r+1)}{\mathrm{\Gamma }(n-\alpha )\mathrm{\Gamma }(r-n+1)}$$ $$\times {\int }_0^t{t}^{n-\alpha -1}{u}^{r-n}{(1-\frac{u}{t})}^{n-\alpha -1}du.$$

Realizando a mudança de variáveis u = ty, tem-se que,

$${\mathbf{D}}_t^{\alpha }[t^r]=\frac{\mathrm{\Gamma }(r+1)t^{r-\alpha }}{\mathrm{\Gamma }(n-\alpha )\mathrm{\Gamma }(r-n+1)}{\int }_0^1{y}^{r-n}{(1-y)}^{n-\alpha -1}\mathit{d}y.$$

Assim,

$${\mathbf{D}}_t^{\alpha }[t^r]=\frac{\mathrm{\Gamma }(r+1)t^{r-\alpha }}{\mathrm{\Gamma }(n-\alpha )\mathrm{\Gamma }(r-n+1)}\beta (r-n+1,n-\alpha ).$$

Temos da Eq. (5) que

$$\beta (r-n+1,n-\alpha )=\frac{\mathrm{\Gamma }(r-n+1)\mathrm{\Gamma }(n-\alpha )}{\mathrm{\Gamma }(r-\alpha +1)},$$

obtêm-se a derivada fracionária segundo Caputo,

(13) $${\mathbf{D}}_t^{\alpha }[t^r]=\frac{\mathrm{\Gamma }(r+1)}{\mathrm{\Gamma }(r-\alpha +1)}{t}^{r-\alpha }.$$

É imediato notar que a derivada fracionária de uma constante, segundo a definição de Caputo, é sempre igual a zero, o que não ocorre no caso da derivada fracionária segundo Riemann-Liouville [¹⁵, ¹⁶], por exemplo. Esse fato é uma das justificativas por se optar pela proposta de Caputo.

Na sequência, apresentamos a transformada de Laplace da derivada fracionária segundo Caputo. A demonstração dessa propriedade será deixada como exercício ao leitor e pode ser encontrada na referência [¹⁵, ¹⁶]. Assim, temos

(14) $$ℒ\{{\mathbf{D}}_t^{\alpha }f(t)\}=s^{\alpha }ℒ\{f(t)\}-\sum _{k=0}^{n-1}{s}^{\alpha -k-1}{f}^{(k)}(0),$$

em que $f^{(k)}(0)$ representa a derivada de ordem inteira k da função $f(t)$ no ponto t=0.

3. Equações Diferenciais Fracionárias

Nesta seção apresentaremos exemplos de equações diferenciais fracionárias no escopo de Caputo que podemos solucionar analiticamente. Não abordaremos teoremas que garantem a existência e unicidade da solução, pois fogem ao escopo desta apresentação, contudo o leitor interessado pode encontrar esse assunto nas referências [¹⁵, ²⁰].

As equações diferenciais fracionárias que serão estudadas neste trabalho possuem a forma geral

(15) $${\mathbf{D}}_t^{\alpha }y(t)=f(t,y(t)),$$ $$y^{(k)}(0)=y_k,$$

em que $k=\mathrm{0,1,2},\mathrm{\dots },\lceil \alpha \rceil -1$, $f(t,y(t))$ é uma função suficientemente regular de t e $y(t)$, e $y^{(k)}(0)$ representa a derivada de ordem inteira k da função $y(t)$ no instante t=0.

Mostraremos a seguir o método da trasformada de Laplace para tratar tais problemas. A viabilidade da solução analítica através desse tipo de estratégia dependerá da forma da função $f(t,y(t))$. Observamos que nas aplicações que serão apresentadas na seção ⁵ apenas o oscilador harmônico permitirá essa análise analítica, pois é a única em que a função f é linear na variável y. No caso da equação logística veremos que a solução analítica através da transformada de Laplace é inviável, considerando que neste caso a $f(t,y)=y(1-y)$ que é não linear.

Nesse sentido, como exemplo ilustrativo da utilização da transformada de Laplace na solução de equações fracionárias, vamos resolver o problema de valor inicial fracionário dado por

(16) $${\mathbf{D}}_t^{\alpha }y=-y,$$ $$y(0)=y_0,$$

em que $0<\alpha <1$. Para esse fim, tomemos a transformada de Laplace da Eq. (16), isto é

$$ℒ\{{\mathbf{D}}_t^{\alpha }y\}=-ℒ\{y\}$$

Denominando $ℒ\{y(t)\}=Y(s)$ e utilizando os resultados das seções anteriores obtemos

$$s^{\alpha }Y(s)-s^{\alpha -1}{y}_0=-Y(s),$$

ou seja

$$Y(s)=\frac{s^{\alpha -1}}{s^{\alpha }+1}{y}_0.$$

Usando a Eq. (11), podemos calcular a transformada inversa da última equação, obtendo

(17) $$y(t)=y_0{E}_{\alpha }(-t^{\alpha }).$$

Ou seja, a solução é dada em termos da função de Mittag-Leffler de um parâmetro. O gráfico apresentado na Figura 1 mostra o comportamento da solução dada na Eq. (17) para diferentes valores de $\alpha$ ($\mathrm{0,5}$; $\mathrm{0,75}$; $\mathrm{0,90}$; $\mathrm{0,95}$ e $\mathrm{1,00}$). Observamos que à medida que aumentamos o valor da ordem da derivada fracionária, o decrescimento da função y se torna mais rápido, lembrando que quando $\alpha =1$ temos que $E_1(-t)=e^{-t}$.

Figura 1
Gráfico da solução do problema de valor inicial Eq. (15) para diferentes valores de $\alpha$ e y0=1.

4. Um Método Numérico para Equações Diferenciais Fracionárias

Diferentemente das derivadas comuns de ordem inteira, que são funcionais locais, as derivadas fracionárias, por definição, são funcionais que possuem uma memória total de estados passados. Isso se deve ao fato das derivadas fracionárias serem definidas a partir de uma convolução entre um núcleo e uma função. Essa observação nos leva a crer que os métodos do tipo Runge-Kutta não serão os mais eficientes para tratar as equações fracionárias, exatamente por levarem em consideração apenas o comportamento local da função. Notamos aqui que métodos dessa forma podem ser facilmente estendidos para os casos de equações diferenciais fracionárias através da série de Taylor generalizada e são encontrados na literatura apresentando bons resultados para casos particulares [¹⁷, ¹⁸, ¹⁹]. Apesar disso, nosso interesse é trabalhar com um método numérico que não tenha as limitações locais dos métodos do tipo Runge-Kutta. Portanto, iremos descrever nessa seção um método em duas etapas do tipo preditor-corretor, que pode ser visto como uma variação do método de Adams-Bashforth, para lidar com problemas do tipo dado pela Eq. (15). O método aqui descrito é devido a Diethelm et al. [¹, ²⁰] e maiores detalhes analíticos sobre o método, como análise de erro e convergência, podem ser encontrados nas referências [¹, ²⁰, ²¹, ²²]. A descrição do método aqui será feita de tal forma que o leitor possa desenvolver métodos numéricos análogos que possam satisfazer suas necessidades.

Aqui queremos tratar do problema de valor inicial dado pela equação diferencial fracionária

(18) $${\mathbf{D}}_t^{\alpha }y(t)=f(t,y(t)),$$

em que $\alpha \in (\mathrm{0,1})$ e $t\in [0,T]$, com a condição inicial

$$y(0)=y_0.$$

Observamos que é possível descrever o método para $\alpha >0$ e a escolha por $\alpha \in (\mathrm{0,1})$ é apenas por motivos didáticos.

Começamos a construção do método numérico definindo o operador integral

$${\mathbf{J}}_t^{\alpha }g(t)=\frac{1}{\mathrm{\Gamma }(\alpha )}{\int }_0^t{(t-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega ,$$

conhecido como integral fracionária de Riemann-Liouville [¹⁶, ²⁰]. Aplicando o operador ${\mathbf{J}}_t^{\alpha }$ em ${\mathbf{D}}_t^{\alpha }y(t)$ temos que

(19) $${\mathbf{J}}_t^{\alpha }{\mathbf{D}}_t^{\alpha }y(t)={\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle {\int }_0^t}{(t-\omega )}^{\alpha -1}{\mathbf{D}}_{\omega }^{\alpha }y(\omega )\mathit{d}\omega$$ $$={\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle {\int }_0^t}{(t-\omega )}^{\alpha -1}$$ $$\times \left[{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\displaystyle {\int }_0^{\omega }}{(\omega -\mu )}^{-\alpha }{y}^{\prime }(\mu )\mathit{d}\mu \right]d\omega$$ $$={\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\displaystyle {\int }_0^t}{(t-\omega )}^{\alpha -1}$$ $$\times \left[{\displaystyle {\int }_0^{\omega }}{(\omega -\mu )}^{-\alpha }{y}^{\prime }(\mu )\mathit{d}\mu \right]d\omega .$$

Fazendo uma mudança de ordem de integração é possível reescrever a Eq. (19) na forma

(20) $${\mathbf{J}}_t^{\alpha }{\mathbf{D}}_t^{\alpha }y(t)=\frac{1}{\mathrm{\Gamma }(\alpha )}\frac{1}{\mathrm{\Gamma }(1-\alpha )}{\int }_0^t{y}^{\prime }(\mu )$$ $$\times \left[{\int }_{\mu }^t{(t-\omega )}^{\alpha -1}{(\omega -\mu )}^{-\alpha }\mathit{d}\omega \right]d\mu .$$

Agora com a mudança de variável $s=\omega -\mu$ podemos verificar que a integral em Eq. (20) é dada por

(21) $${\mathbf{J}}_t^{\alpha }{\mathbf{D}}_t^{\alpha }y(t)$$ $$=\frac{1}{\mathrm{\Gamma }(\alpha )}\frac{1}{\mathrm{\Gamma }(1-\alpha )}{\int }_0^t{y}^{\prime }(\mu )$$ $$\times \left[{\int }_0^{t-\mu }{s}^{-\alpha }{(t-\mu -s)}^{\alpha -1}\mathit{d}s\right]d\mu$$ $$=\frac{1}{\mathrm{\Gamma }(\alpha )}\frac{1}{\mathrm{\Gamma }(1-\alpha )}{\int }_0^t{y}^{\prime }(\mu )$$ $$\times \left[{\int }_0^{t-\mu }{s}^{-\alpha }{(t-\mu )}^{\alpha -1}{\left(1-\frac{s}{t-\mu }\right)}^{\alpha -1}\mathit{d}s\right]d\mu .$$

Novamente realizando uma mudança de variável fazendo $r=t/(t-\mu )$ temos que a integral em Eq. (21) será reescrita na forma

$${\mathbf{J}}_t^{\alpha }{\mathbf{D}}_t^{\alpha }y(t)=\frac{1}{\mathrm{\Gamma }(\alpha )}\frac{1}{\mathrm{\Gamma }(1-\alpha )}{\int }_0^t{y}^{\prime }(\mu )$$ $$\times \left[{\int }_0^1{r}^{(1-\alpha )-1}{\left(1-r\right)}^{\alpha -1}\mathit{d}r\right]d\mu .$$

Notamos, pela identidade dada na Eq. (5), que

$$\beta (\alpha \mathrm{,1}-\alpha )=\mathrm{\Gamma }(\alpha )\mathrm{\Gamma }(1-\alpha )={\int }_0^1{r}^{(1-\alpha )-1}{\left(1-r\right)}^{\alpha -1}\mathit{d}r$$

e finalmente temos o efeito da aplicação do operador ${\mathbf{J}}_t^{\alpha }$ na derivada fracionária ${\mathbf{D}}_t^{\alpha }y(t)$:

$${\mathbf{J}}_t^{\alpha }{\mathbf{D}}_t^{\alpha }y(t)={\int }_0^t{y}^{\prime }(\mu )\mathit{d}\mu$$ $$=y(t)-y(0).$$

Aplicando o operador ${\mathbf{J}}_t^{\alpha }$ na equação diferencial fracionária dada pela Eq. (18) encontramos a seguinte equação integral que será a base para o método aqui descrito

$$y(t)=y(0)+\frac{1}{\mathrm{\Gamma }(\alpha )}{\int }_0^t{(t-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega ,$$

em que $g(\omega )=f(\omega ,y(\omega ))$.

Considere $\{0=t_0<t_1<t_2<\mathrm{\dots }<t_n<t_{n+1}=T\}$ uma discretização equiespaçada do domínio $[0,T]$ com $t_k-t_{k-1}=h$, para todo $k\in \{\mathrm{1,2,3},\mathrm{\dots },n+1\}$. Assim, no ponto $t=t_{k+1}$ temos

(22) $$y_{k+1}=y_0+{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle {\int }_0^{t_{k+1}}}{(t_{k+1}-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega$$ $$=y_0+{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle \sum _{j=0}^k}{\displaystyle {\int }_{t_j}^{t_{j+1}}}{(t_{k+1}-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega ,$$

em que $y(t_k)=y_k,$ para todo $k\in \{\mathrm{0,1,2},\mathrm{\dots },n+1\}$. O que faremos agora é tentar aproximar as integrais que aparecem na Eq. (22) usando a regra dos trapézios, considerando o termo ${(t_{k+1}-\cdot )}^{\alpha -1}$ como uma função peso. Logo, sendo $g_{j+1}(\omega )$ a interpolação linear da função $g(\omega )$ no intervalo $[t_j,t_{j+1}]$ temos que

$$g_{j+1}(\omega )=\frac{\omega -t_j}{h}g(t_{j+1})-\frac{\omega -t_{j+1}}{h}g(t_j).$$

Assim, as integrais no lado direito da Eq. (22) podem ser aproximadas por

$$\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{(t_{k+1}-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega$$ $$\approx \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{(t_{k+1}-\omega )}^{\alpha -1}{g}_{j+1}(\omega )\mathit{d}\omega$$ $$=\sum _{j=0}^{k+1}{\phi }_{j,k+1}g(t_j),$$

em que os coeficientes ${\phi }_{j}k+1$ são calculadas através das seguintes integrais

(23) $${\phi }_{0,k+1}=-{\displaystyle \frac{1}{h}}{\displaystyle {\int }_{t_0}^{t_1}}{(t_{k+1}-\omega )}^{\alpha -1}\left(\omega -t_1\right)\mathit{d}\omega ,$$

(24) $${\phi }_{j,k+1}={\displaystyle \frac{1}{h}}{\displaystyle {\int }_{t_{j-1}}^{t_j}}{(t_{k+1}-\omega )}^{\alpha -1}\left(\omega -t_{j-1}\right)\mathit{d}\omega$$ $$-{\displaystyle \frac{1}{h}}{\displaystyle {\int }_{t_j}^{t_{j+1}}}{(t_{k+1}-\omega )}^{\alpha -1}\left(\omega -t_{j+1}\right)\mathit{d}\omega ,$$ $$j\in \{1,\mathrm{\dots },k\}$$

(25) $${\phi }_{k+1,k+1}={\displaystyle \frac{1}{h}}{\displaystyle {\int }_{t_k}^{t_{k+1}}}{(t_{k+1}-\omega )}^{\alpha -1}\left(\omega -t_k\right)\mathit{d}\omega .$$

Observamos que as integrais dadas pelas Eqs. (23)–(25) são imediatas e calculando-as obtemos os coeficientes ${\phi }_{j}k+1$ desejados

(26) $${\phi }_{0,k+1}=\frac{h^{\alpha }}{\alpha (\alpha +1)}\left[k^{\alpha +1}-(k-\alpha ){(k+1)}^{\alpha }\right],$$

(27) $${\phi }_{j,k+1}=\frac{h^{\alpha }}{\alpha (\alpha +1)}[{(k-j+2)}^{\alpha +1}+{(k-j)}^{\alpha +1}$$ $$-2{(k-j+1)}^{\alpha +1}],j\in \{1,\mathrm{\dots },k\},$$

(28) $${\phi }_{k+1,k+1}=\frac{h^{\alpha }}{\alpha (\alpha +1)}.$$

Finalmente, podemos escrever que a solução aproximada no instante de tempo $t=t_{k+1}$ é dada por

(29) $$y_{k+1}=y_0+{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\displaystyle \sum _{j=0}^k}{\phi }_{j,k+1}f(t_j,y_j)$$ $$+{\displaystyle \frac{1}{\mathrm{\Gamma }(\alpha )}}{\phi }_{k+1,k+1}f(t_{k+1},y_{k+1}^p).$$

Note que na fórmula encontrada na Eq. (29), para a aproximação da solução da equação diferencial, o valor da função $y(t)$ que queremos calcular no instante $t_{k+1}$, $y_{k+1}$, aparece naturalmente também no lado direito da expressão, como argumento da função f. Como esse valor não é conhecido precisamos estimá-lo e ao valor desta estimativa denominaremos $y_{k+1}^p$. Para tal utilizaremos um procedimento análogo ao já descrito anteriormente modificando apenas a quadratura realizada nas integrais da Eq. (22), ao invés de as aproximarmos utilizando a regra dos trapézios iremos considerar que a função dentro de cada subintervalo $[t_j,t_{j+1}]$ é constante e igual a $g(t_j)$. Isto é,

$$\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{(t_{k+1}-\omega )}^{\alpha -1}g(\omega )\mathit{d}\omega$$ $$\approx \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{(t_{k+1}-\omega )}^{\alpha -1}{q}_{j+1}(\omega )\mathit{d}\omega$$ $$=\sum _{j=0}^k{\varrho }_{j,k+1}q(t_j),$$

com $q_{j+1}(\omega )=g(t_j),$ para todo $\omega \in [t_j,t_{j+1}]$, sendo os coeficientes ${\varrho }_{j,k+1}$ dados por

(30) $${\varrho }_{j,k+1}={\int }_{t_j}^{t_{j+1}}{(t_{k+1}-\omega )}^{\alpha -1}\mathit{d}\omega$$ $$=\frac{h^{\alpha }}{\alpha }\left[{(k-j+1)}^{\alpha }-{(k-j)}^{\alpha }\right].$$

Finalmente, temos a partir da Eq. (22) que a estimativa para $y_{k+1}^p$ pode ser escrita como

(31) $$y_{k+1}^p=y_0+\frac{1}{\mathrm{\Gamma }(\alpha )}\sum _{j=0}^k{\varrho }_{j,k+1}f(t_j,y_j).$$

Portanto, a solução numérica a ser obtida pelo método de Adam-Bashforth fracionário é construída a partir da predição de $y_{k+1}^p$, utilizando a Eq. (31), seguida da correção dada pela Eq. (29) e que fornece a aproximação desejada $y_{k+1}$ correspondente à $y(t_{k+1})$.

A fim de verificarmos a qualidade do método aqui descrito, temos no gráfico apresentado na Figura 2 uma comparação entre a solução analítica e numérica para um problema de valor inicial. Aqui foi determinada a solução do problema de valor inicial dado pela equação ${\mathbf{D}}_t^{1/2}=-y$, no intervalo $[\mathrm{0,2}]$, com condição inicial $y(0)=1$. Conforme mostramos na Eq. (17) a solução analítica deste problema é dada pela função de Mittag-Leffler $y(t)=E_{1/2}(-t^{1/2})$. No caso da solução numérica observamos que foi utilizado o algoritmo de Adam-Bashforth fracionário com passo $h=5\times {10}^{-2}$. É imediato verificar qualitativamente a adequação da solução numérica em relação a solução analítica.

Finalmente, se faz importante observar que para problemas de valor inicial para $\alpha >1$ na forma

$${\mathbf{D}}_t^{\alpha }y(t)=f(t,y(t)),$$ $$y^{(m)}(0)=y_0^{(m)},m=\mathrm{0,1},\mathrm{\dots },\lceil \alpha \rceil -1$$

é imediato mostrar que as Eqs. (29)–(31) são ligeiramente modificadas e assumem a forma

(32) $$y_{k+1}=\sum _{m=0}^{\lceil \alpha \rceil -1}\frac{t_{k+1}^m}{m!}{y}_0^{(m)}+\frac{1}{\mathrm{\Gamma }(\alpha )}\sum _{j=0}^k{\phi }_{j,k+1}f(t_j,y_j)$$ $$+\frac{1}{\mathrm{\Gamma }(\alpha )}{\phi }_{k+1,k+1}f(t_{k+1},y_{k+1}^p)$$

(33) $$y_{k+1}^p=\sum _{m=0}^{\lceil \alpha \rceil -1}\frac{t_{k+1}^m}{m!}{y}_0^{(m)}+\frac{1}{\mathrm{\Gamma }(\alpha )}\sum _{j=0}^k{\varrho }_{j,k+1}f(t_j,y_j).$$

Este formato será importante na aplicação do oscilador harmônico fracionário.

Figura 2
Solução analítica e numérica da equação ${\mathbf{D}}_t^{1/2}=-y$ no intervalo $[\mathrm{0,2}]$ com condição inicial $y(0)=1$. Para este problema a solução analítica é dada pela função de Mittag-Leffler $y(t)=E_{1/2}(-t^{1/2})$. Além disso, para a solução numérica foi utilizado o passo $h=5\times {10}^{-2}$.

5. Aplicações

Nesta seção aplicaremos o método descrito na seção anterior para analisar brevemente três modelos baseados em equações diferenciais fracionárias. O primeiro deles é o oscilador harmônico fracionário submetido a uma força externa, outro é a versão fracionária da equação logística e, por fim, trataremos de um modelo de viscoelasticidade para tecidos moles. Esperamos que o tratamento aqui apresentado contribua de forma significativa nas pesquisas sobre cálculo fracionário e sirva de referência literária introdutória na área.

5.1. O oscilador harmônico fracionário submetido a uma força externa senoidal

O oscilador harmônico é um sistema físico constituído por um objeto de massa m preso a uma mola de constante k. A equação diferencial que governa o movimento do oscilador, isto é, a equação que fornece a posição do oscilador em função do tempo é dada por

(34) $$\frac{d^{2}x(t)}{d{t}^2}+{\omega }_o^{2}x(t)=0,$$

em que ${\omega }_o^2=\frac{k}{m}$. A Eq. (34) é usada para analisar o oscilador livre. Se por acaso o oscilador estiver submetido a alguma força externa de módulo F_ext, podemos escrever

(35) $$\frac{d^{2}x(t)}{d{t}^2}+{\omega }_o^{2}x(t)=F_{ext}.$$

É relevante aqui recordarmos que o oscilador harmônico é um modelo que descreve bem muitos sistemas físicos reais, como por exemplo sistemas que oscilam em torno de uma posição de equilíbrio. Ademais, se considerarmos um oscilador com dissipação, isto é, se as forças de atrito forem consideradas, obtemos os casos especiais do oscilador amortecido, superamortecido e subamortecido, cada qual com a sua aplicação prática. Um fato bastante peculiar no caso do oscilador harmônico fracionário é a possibilidade de se descrever dissipação sem a necessidade de introduzir o termo associado à força de atrito na equação do movimento; bastando para isso variar a ordem da derivada fracionária. Além disso, podemos generalizar o oscilador para um tratamento quântico, submetendo a equação de Schrödinger ao potencial do oscilador harmônico, o que nos permite descrever por exemplo o comportamento de moléculas diatômicas com bastante acurácia. Porém, o problema do oscilador harmônico no âmbito da mecânica quântica ainda é um problema muito pouco explorado.

O correspondente fracionário do modelo físico dado na Eq. (35) quando submetido a uma força externa senoidal é dado por [³, ³⁷]

(36) $${\mathbf{D}}_t^{\alpha }x(t)+{\omega }_o^{\alpha }x(t)=F_o\sin (\omega t),$$

em que $\alpha \in (\mathrm{1,2})$ e F_o é a amplitude da força externa. Observe que quando $1<\alpha <2$ duas condições iniciais são necessárias para a solução do problema, a saber: $x(0)$ e $x^{\prime }(0)$. Observa-se que uma interpretação física para o oscilador harmônico fracionário foi dada por Stanislavsky [²³]. Ele sugeriu que um oscilador harmônico fracionário deve ser interpretado como uma média de osciladores harmônicos comuns (de derivada ordinária) governado por uma seta do tempo estocástica. Além disso, do ponto de vista matemático os osciladores harmônicos fracionários também são úteis para determinar a solução da equação de Schrödinger fracionária livre [⁷, ²⁴]. Ressalta-se ainda que aspectos fisicamente significativos do sistema dinâmico proporcionado pela equação do oscilador harmônico fracionário, como energia total e resposta dinâmica do sistema, também já foram estudados na literatura [²⁵, ²⁶, ²⁷].

A solução da Eq. (36) para o caso livre, isto é, F₀ = 0, pode ser obtida de maneira analítica diretamente a partir da transformada de Laplace, considerando as Eq. (11) e Eq. (14), e é dada por

(37) $$x(t)=x(0)E_{\alpha \mathrm{,1}}(-{\omega }_0^{\alpha }{t}^{\alpha })+x^{\prime }(0)t{E}_{\alpha \mathrm{,2}}(-{\omega }_0^{\alpha }{t}^{\alpha }).$$

No entanto, o mesmo não se pode dizer sobre a solução analítica para o caso em que $F_0\ne 0$. Ainda é possível encontrar solução analítica para casos particulares da Eq. (36) utilizando a transformada de Laplace e o teorema dos resíduos, mas a solução está longe de ser encontrada trivialmente. Por exemplo, para $\alpha =3/2$, ${\omega }_0=\omega =1$, F⁰ = 1 e $x(0)=x^{\prime }(0)=1$ é possível mostrar que a solução analítica da Eq. (36) é

(38) $$x(t)=E_{\frac{3}{2}\mathrm{,1}}(-t^{3/2})+t{E}_{\frac{3}{2}\mathrm{,2}}(-t^{3/2})+{\displaystyle \frac{4}{3}}{e}^{-t/2}\cos \left({\displaystyle \frac{\sqrt{3}}{2}}t\right)$$ $$+{\displaystyle \frac{\text{sen}(t)}{2}}-{\displaystyle \frac{\cos (t)}{2\sqrt{2}-2}}.$$

Nesse sentido, trataremos aqui a Eq. (36) utilizando o método numérico desenvolvido na seção anterior. Nosso objetivo aqui é duplo: no primeiro momento estudaremos o que ocorre no modelo fracionário quando $\omega \ne {\omega }_o$, e em seguida, veremos como o sistema fracionário se comporta no caso em que a frequência natural coincide com a da força externa, isto é, $\omega ={\omega }_o$.

Primeiramente, vamos avaliar a qualidade da aproximação numérica a comparando à solução analítica dada pela Eq. (38). Na Figura 3 podemos ver tal comparação e qualitativamente é possível observar que a solução obtida com o método numérico descrito se aproxima satisfatoriamente da solução analítica.

Figura 3
Solução analítica e numérica da equação ${\mathbf{D}}_t^{3/2}=\mathrm{sen}(t)-y$ no intervalo $[\mathrm{0,100}]$ com condições iniciais $y(0)=1$ e $y^{\prime }(0)=1$. Para este problema a solução analítica é dada pela Eq. (38). Além disso, para a solução numérica foi utilizado o passo $h={10}^{-1}$.

Na Figura 4 podemos ver o caso do oscilador harmônico para ${\omega }_0={(\mathrm{0,5})}^{2/3}$ e $\omega =1$ para diferentes valores de $\alpha$ ($\mathrm{1,2}$; $\mathrm{1,4}$; $\mathrm{1,6}$ e $\mathrm{1,8}$). Percebemos que quanto menor é a ordem da derivada fracionária, menor é a amplitude da oscilação. E ainda, percebemos também que neste caso há uma espécie de amortecimento da oscilação, mesmo não tendo considerado ações de forças dissipativas.

Figura 4
Solução numérica do oscilador harmônico fracionário ${\mathbf{D}}_t^{\alpha }=\mathrm{sen}(t)-\frac{1}{2}y$ no intervalo $[\mathrm{0,100}]$ com condições iniciais $y(0)=1$ e $y^{\prime }(0)=-1$ para diferentes valores de $\alpha \in (\mathrm{1,2})$. Além disso, para a solução numérica foi utilizado o passo $h={10}^{-1}$.

Finalmente, na Figura 5 podemos ver o caso do oscilador harmônico para ${\omega }_0=\omega =1$ para diferentes valores de $\alpha$ ($\mathrm{1,2}$; $\mathrm{1,4}$; $\mathrm{1,6}$ e $\mathrm{1,8}$). Conforme já conhecido da literatura, quando as frequência da fonte externa é igual à frequência natural de oscilação do sistema, no cálculo de ordem inteira, temos o fenômeno da ressonância, ou seja, a amplitude vai crescendo indefinidamente. Contudo, notamos no gráfico e também no exemplo dado pela equação (38), que enquanto a ordem da derivada fracionária diminui, a oscilação se estabiliza, isto é, além da amplitude do movimento oscilatório ser menor, a ressonância deixa de acontecer. No caso de oscilações forçadas de um sistema mecânico, infere-se deste resultado que mesmo uma força externa com mesma frequência não levará o sistema a se romper.

Figura 5
Solução numérica do oscilador harmônico fracionário ${\mathbf{D}}_t^{\alpha }=\mathrm{sen}(t)-y$ no intervalo $[\mathrm{0,100}]$ com condições iniciais $y(0)=1$ e $y^{\prime }(0)=-1$ para diferentes valores de $\alpha \in (\mathrm{1,2})$. Além disso, para a solução numérica foi utilizado o passo $h={10}^{-1}$.

5.2. Equação logística fracionária

A equação logística

(39) $$\frac{dy(t)}{dt}=ky(t)(1-y(t)),$$

em que k é uma constante, foi proposta originalmente por Verhulst em 1838 [²⁸]. Em seu trabalho, Verhulst tentou aprimorar o modelo de Malthus para a estatística populacional [²⁹]. Desde então, a equação logística passou a ser aplicada em diversos contextos, desde modelos biológicos de proliferação de bactérias até predição de sistemas econômicos [³⁰, ³¹].

Uma hipótese básica relacionada à Eq. (39) é que, para o processo de variação da população, o tempo é uma quantidade passiva homogênea que indexa os eventos de mudança da quantidade por meio do tempo regularmente medido por um relógio. No entanto, este não é necessariamente o caso para fenômenos complexos e, assim, a derivada temporal é continuada analiticamente para uma derivada fracionária para produzir a equação logística de ordem fracionária [³², ³⁰]

(40) $${\mathbf{D}}_t^{\alpha }y(t)=k^{\alpha }y(t)(1-y(t)),$$

em que k é uma constante, $0<\alpha <1$ e com uma condição inicial $y(0)=y_0$. Cabe ressaltar que devido aos efeitos de memória e também a não-localidade das derivadas de Caputo, diversos modelos populacionais da literatura, e não apenas a equação logística, recentemente foram estendidos para suas respectivas versões fracionárias com o intuito de avaliar melhor a evolução de doenças [³³, ³⁴, ³⁵].

Figura 6
Solução numérica da equação logística fracionária ${\mathbf{D}}_t^{\alpha }=y(1-y)$ no intervalo $[\mathrm{0,10}]$ para diferentes valores de $\alpha \in (\mathrm{0,1})$. (a) Solução com condição inicial $y_0=1/2$; (b) solução com condição inicial $y_0=3/2$. Além disso, para a solução numérica foi utilizado o passo (a) $h={10}^{-2}$; e (b) $h={10}^{-3}$.

Note que diferentemente do oscilador harmônico a Eq. (40) é não-linear e portanto os métodos envolvendo a transformada de Laplace dada pela Eq. (14) não podem ser aplicados. A resolução desta equação pode aparecer na literatura de algumas formas distintas: o tratamento de uma versão linearizada [³¹], solução analítica através de expansão em algum tipo de série [³²], solução analítica através de alguma hipótese específica sobre a forma da solução [³⁶] e abordagem via métodos numéricos [³⁰]. A nossa abordagem aqui será realizada com o uso do método construído na seção precedente.

Na Figura 6(a) podemos ver as soluções da equação logística para diferentes valores de $\alpha$ ($\mathrm{0,2}$; $\mathrm{0,4}$; $\mathrm{0,6}$; $\mathrm{0,8}$ e $\mathrm{1,0}$) com a condição inicial $y_0=1/2$ e $h={10}^{-2}$. Vemos assim que quanto menor a ordem da derivada fracionária, a solução da equação logística fracionária cresce de forma mais suave e se estabiliza num valor mais baixo, isto é, o valor suporte é menor. Se tomarmos como exemplo o crescimento de uma população, se o problema puder ser modelado via cálculo fracionário, a população que cresce mais devagar e se estabiliza num valor menor é modelada por uma derivada de ordem fracionária mais baixa.

Por outro lado, na Figura 6(b) vemos as soluções da equação logística para diferentes valores de $\alpha$ ($\mathrm{0,2}$; $\mathrm{0,4}$; $\mathrm{0,6}$; $\mathrm{0,8}$ e $\mathrm{1,0}$) com a condição inicial $y_0=3/2$ e $h={10}^{-3}$. Note que neste caso o passo é menor que o exemplo anterior para garantir uma solução suave para o caso $\alpha =\mathrm{0,2}$. Novamente podemos ver que há uma estabilização mais suave da solução quando a ordem da derivada fracionária é menor.

5.3. Modelo de viscoelasticidade de ordem fracionária para tecidos moles

Os tecidos biológicos moles podem ser considerados como materiais viscoelásticos, isto é, apresentam propriedades elásticas e viscosas ao mesmo tempo [³⁸].

A lei elástica de Fung [³⁹] é o modelo padrão-ouro na literatura para modelar a resposta elástica de tecidos moles [³⁸, ¹¹, ¹²]. Ela descreve o comportamento não-linear da tensão-deformação de um tecido biológico mole e foi deduzida a partir de observações experimentais.

A lei de Fung é descrita pela equação diferencial ordinária

(41) $$\frac{ds(\lambda )}{d\lambda }=E+\beta s(\lambda ),s(1)=0,$$

em que s é a tensão Lagrangiana, $\lambda$ o alongamento, E é uma constante com unidades de tensão e $\beta$ é uma constante adimensional. Ambas as constantes são características do material envolvido.

De forma não-trivial pode-se estender a lei de Fung para o modelo [¹¹]

(42) $$\left(1+\tau \frac{d}{dt}\right)s(t)=\frac{E}{\beta }\left(e^{\beta ϵ(t)}-1\right)+\rho (E+\beta s(t))\frac{d}{dt}ϵ(t)$$

que descreve as características viscoelásticas dos tecidos moles. Aqui os parâmetros materiais são constantes não-negativas tais que E tem unidades de tensão, $\beta$ é adimensional, $\rho$ e $\tau$ tem unidades de tempo com $\rho >\tau$. Além disso, a função $ϵ(t)$ representa a deformação e é supostamente conhecida. No entanto, este modelo apresentado peca ao modelar com precisão dados experimentais [¹¹].

Já foi mostrado que as propriedades viscoelásticas desses tecidos são melhores tratadas com modelos viscoelásticos de ordem fracionária [⁴⁰, ⁴¹]. Desta forma, o modelo dado pela Eq. (42) é continuado analiticamente para a equação diferencial fracionária [¹¹]

(43) $$\left(1+{\tau }^{\alpha }{\mathbf{D}}_t^{\alpha }\right)s(t)=\frac{E}{\beta }\left(e^{\beta ϵ(t)}-1\right)+{\rho }^{\alpha }(E+\beta s(t)){\mathbf{D}}_t^{\alpha }ϵ(t),$$

em que $0<\alpha <1$.

Aqui vamos replicar o experimento do teste de tensão da referência [¹¹]. Com isso, iremos simular o comportamento do colágeno para uma função de deformação $ϵ(t)$ conhecida. Desta forma, a Eq. (43) será resolvida numericamente com os seguintes parâmetros: $\alpha =\mathrm{0,33}$; $\beta =\mathrm{12,4}$; E=1 MPa, $\tau =10\mu$s e $\rho =10$ ms. O teste de tensão será feito utilizando uma deformação linear $ϵ(t)=kt$, com taxa de deformação k constante, para um material em estado virgem, isto é, $s(0)=0$. Além disso, serão considerados três valores distintos de k, a saber: $k=\mathrm{0,01}$, $k=\mathrm{0,01}$ e $k=\mathrm{0,001}$. A simulação será feita no intervalo de tempo $0\le t\le t_{\max }$, com t_max escolhido tal que $ϵ(t_{\max })=\mathrm{0,15}$. Finalmente, para a função de deformação escolhida temos, pela Eq. (13), que ${\mathbf{D}}_t^{\alpha }ϵ(t)=\frac{k{t}^{1-\alpha }}{\mathrm{\Gamma }(2-\alpha )}$. Com isso, o problema a ser resolvido numericamente é dado pela equação fracionária linear

(44) $${\mathbf{D}}_t^{\alpha }s(t)={\tau }^{-\alpha }\left({\rho }^{\alpha }\beta {\displaystyle \frac{k{t}^{1-\alpha }}{\mathrm{\Gamma }(2-\alpha )}}-1\right)s(t)$$ $$+{\tau }^{-\alpha }E\left({\displaystyle \frac{1}{\beta }}\left(e^{\beta kt}-1\right)+{\rho }^{\alpha }{\displaystyle \frac{k{t}^{1-\alpha }}{\mathrm{\Gamma }(2-\alpha )}}\right).$$

Na Figura 7 podemos ver a solução para o teste de tensão com gráfico apresentado na forma tensão-deformação. Os resultados apresentados são compatíveis com resultados experimentais e correspondem aos mesmos apresentados por Freed et al. [¹¹]. Novamente, percebe-se a adequação do método numérico ao problema de viscoelasticidade do colágeno e a sua importância para a solução de problemas teóricos e reais.

Figura 7
Resultado do teste de tensão para a simulação do colágeno modelado pela Eq. (44). Aqui foram utilizados os seguintes dados: deformação linear $ϵ(t)=kt$, com $k=\mathrm{0,1}$, $k=\mathrm{0,01}$ e $k=\mathrm{0,001}$. Cada uma das simulações foi realizada no intervalo de tempo $[0,t_{\max }]$, de tal forma que para cada valor de k a igualdade $ϵ(t_{\max })=\mathrm{0,15}$ deve ser respeitada.

6. Considerações Finais e Perspectivas

Neste trabalho, apresentam-se os conceitos das equações diferenciais fracionárias no escopo de Caputo enfatizando sua resposta analítica que é deduzida pelo método da Transformada de Laplace cuja solução geral é dada em termos da função de Mittag-Leffer de um parâmetro.

Além disso, enfatizou-se um método do tipo preditor-corretor para a resolução numérica destas equações diferenciais fracionárias, na qual é uma variação do método de Adam-Bashforth, que caracterizou-se pela sua eficácia em relação aos métodos do tipo Runge-Kutta usualmente abordados na literatura para a solução numérica destas equações.

As três aplicações expostas: o oscilador harmônico fracionário submetido a uma força externa senoidal, a equação logística e o modelo de viscoelasticidade, exemplificam o funcionamento do método numérico descrito, na qual os resultados numéricos obtidos mostraram valores com convergência rápida do valor exato e com pequenos erros.

Espera-se com este trabalho contribuia para uma melhor compreensão das equações diferenciais fracionárias no escopo de Caputo e uma reflexão sobre sua importância aplicação nos fenômenos físicos, tanto de teoria (que traz uma novo conceito teórico) quanto na solução (que traz uma abordagem eficaz para seu desenvolvimento numérico).

Referências

• ^[1]
  K. Diethelm, N.J. Ford e A.D. Freed, Nonlinear Dynamics 29, 3 (2002).
• ^[2]
  B. Ross, Historia Mathematica 4, 75 (1977).
• ^[3]
  R.F. Camargo e E.C. Oliveira, Cálculo fracionário (Livraria da Física, São Paulo, 2015).
• ^[4]
  M. Caputo, Geophysical Journal of the Royal Astronomical Society 13, 529 (1967).
• ^[5]
  K.B. Oldhan e J. Spanier, The Fractional Calculus: Theory and Aplications, Differentiation and Integration to Arbitrary Order (Academic, New York, 1974).
• ^[6]
  N. Varalta, A.V. Gomes e R.F. Camargo, TEMA 15, 211 (2014).
• ^[7]
  R. Herrmann, Fractional Calculus: An Introduction to Physicists (World Scientific, Singapura (2011).
• ^[8]
  P.J. Torvik e R.L. Bagley, Journal of Applied Mechanics 51, 94 (1984).
• ^[9]
  W. Labecca, O. Guimarães e J.R.C. Piqueira, Mathematical Problems in Engineering, 591715 (2015).
• ^[10]
  J. Vaz Jr e E.C. Oliveira, Mathematics in Engineering 4, 1 (2022).
• ^[11]
  A.D. Freed e K. Diethelm, Fractional Calculus and Applied Analysis 10, 219 (2007).
• ^[12]
  T.W.J. Almeida, Sistema para análise viscoelástica de tecidos moles por ondas de cisalhamento usando excitação magnética e medida ultrassônica Tese de Doutorado, Universidade de São Paulo, Ribeirão Preto (2015).
• ^[13]
  J.P.S. Ramírez, Função Gama, disponível em: https://www.ime.unicamp.br/~ftorres/ENSINO/MONOGRAFIAS/JP_VC2_2015.pdf (2015).
  » https://www.ime.unicamp.br/~ftorres/ENSINO/MONOGRAFIAS/JP_VC2_2015.pdf
• ^[14]
  D.S. Oliveira, Derivada Fracionária e as Funções de Mittag-Leffler Dissertação de Mestrado, Universidade Estadual de Campinas, Campinas (2014).
• ^[15]
  A.A. Kilbas, H.M. Srivastava e J.H. Trujillo, Theory and Application of Fractional Differential Equations (Elsevier, Reino Unido, 2006).
• ^[16]
  K.M. Owolabi, A. Atangana, Numerical Methods for Fractional Differentiation (Springer, Singapura, 2019).
• ^[17]
  Z.M. Odibat e N.T. Shawagfeh, Applied Mathematics and Computation 186, 286 (2007).
• ^[18]
  S. ÖLrekçi, Advances in Mathematical Physics, 507970 (2015).
• ^[19]
  M.S. Arshad, D. Baleanu, M.B. Riaz e M. Abbas, Discrete Dynamics in Nature and Society, 1020472 (2020).
• ^[20]
  K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type (Springer, Singapura, 2010).
• ^[21]
  K. Diethelm, N.J. Ford e A.D. Freed, Numerical Algorithms 36, 31 (2004).
• ^[22]
  K. Diethelm, N.J. Ford, A.D. Freed e Y. Luchko, Computer Methods in Applied Mechanics and Engineering, 194, 743 (2005).
• ^[23]
  A.A. Stanislavsky, Physical Review E 70, 051103 (2004).
• ^[24]
  X. Guo e M. Xu, Journal of Mathematical Physics 47, 082104 (2006).
• ^[25]
  B.N.N. Achar, J.W. Hanneken, T. Enck e T. Clarke, Physica A 309, 275 (2002).
• ^[26]
  B.N.N. Achar, J.W. Hanneken, T. Enck e T. Clarke, Physica A 297, 361 (2001).
• ^[27]
  A. Al-rabtah, V.S. Ertürk e S. Momani, Computers & Mathematics with Applications 59, 1356 (2010).
• ^[28]
  P.F. Verhulst, Correspondance mathématique et physique 59, 113 (1838).
• ^[29]
  R.C. Bassanezzi, Ciência e Natura 36, 97 (2014).
• ^[30]
  A.M.A. El-Sayed, A.E.M. El-Mesiry e H.A.A. El-Saka, Applied Mathematics Letters 20, 817 (2007).
• ^[31]
  R.F. Camargo, M.M. Theodoro, Revista Eletrônica Paulista de Matemática 17, 61 (2020).
• ^[32]
  B.J. West, Physica A: Statistical Mechanics and its Applications 429, 103 (2015).
• ^[33]
  I. Ahmed, I.A. Baba, A. Yusuf, P. Kumam, W. Kumam, Advances Differences Equations, 394, 1 (2020).
• ^[34]
  M.A.A. Oud, A. Ali, H. Alrabaiah, S. Ullah, M.A. Khan, A.A. Islam, Advances Differences Equations, 106, 1 (2021).
• ^[35]
  F. Alcántara-López, C. Fuentes, C. Chávez, F. Brambila-Paz, A. Quevedo, Mathematics 9 1915 (2021).
• ^[36]
  V.E. Tarasov, Mathematics 8, 2231 (2020).
• ^[37]
  G.S. Teodoro, D.S. Olievira, C. Oliveira, Rev. Bras. de Ensino de Física 40, e2307 (2018).
• ^[38]
  Y.C. Fung, Biomechanics: mechanical properties of living tissues (Springer Verlag, New York, 1993).
• ^[39]
  Y.C. Fung, American Journal of Physiology 28, 1532 (1967).
• ^[40]
  T.C. Doehring, A.D. Freed, E.O. Carew, I. Vesely, Journal of Biomechanical Engineering 127, 700 (2005).
• ^[41]
  A.D. Freed, K. Diethelm, Biomechanics and Modeling in Mechanobiology 5, 203 (2006).

Datas de Publicação

• Publicação nesta coleção
  17 Jan 2022
• Data do Fascículo
  2022

Histórico

• Recebido
  09 Dez 2021
• Aceito
  15 Dez 2021