Construindo transformadas finitas usando a Teoria de Sturm--Liouville

Setting up finite transforms using Sturm--Liouville Theory

Resumo

Problemas de valores inicial e de contorno são muito comuns na Física, Matemática e Engenharia. Eles podem modelar diversos tipos de problemas relacionados a difusão de calor e a vibração de membranas, por exemplo. Quando se deseja encontrar a solução analítica desses problemas podemos encontrar dificuldades extras quando as equações e também as condições de contorno que descrevem os fenômenos são não-homogêneas. Desta forma, neste trabalho apresentamos uma técnica de solução de problemas de valores iniciais e de contorno por meio de transformações integrais. O diferencial da apresentação está na construção da transformada integral apropriada à solução do problema. Essas transformadas são conhecidas como transformadas finitas e neste caso elas estão relacionadas a um problema de Sturm–Liouville associado com o operador diferencial ligado à equação diferencial. Como exemplo do desenvolvimento e aplicação da ferramenta, resolvemos dois problemas de difusão de calor em coordenadas espaciais distintas. A apresentação do trabalho segue de forma pedagógica e autocontida. Sendo assim, esperamos que o leitor compreenda a técnica e possa utilizá-la na resolução de outros problemas envolvendo equações diferencias parciais.

Palavras-chave:
Equação do Calor; Problema de Sturm--Liouville; Transformada Integral; Transformada Finita

Abstract

Initial and boundary value problems are very common in physics, mathematics, and engineering. They can model many types of problems related to heat diffusion and membrane vibration, for example. When the analytical solution of these problems is needed one may find extra difficulties when the boundary conditions that describe the phenomena are nonhomogeneous. In this way, in this work is presented a technique for solving initial and boundary problems by means of integral transformations. The differential of the presentation is in the construction of the integral transform appropriate to the solution of the problem. These transformations are known as finite transforms and in this case they are related to a Sturm--Liouville problem associated with the differential operator connected to the differential equation. As an example, we solve two problems of heat diffusion in different spatial coordinates systems. The presentation of the work follows in pedagogically and self-contained fashion. Therefore, we expect the reader to understand the technique and can use it in solving other problems involving partial differential equations.

Keywords:
Heat Equation; Sturm--Liouville Problem; Integral Transform; Finite Transform

1. Introdução

De um ponto de vista histórico, podemos afirmar que as equações diferenciais surgiram a partir do século XVII, juntamente com o cálculo diferencial e integral, no contexto das aplicações da matemática em física e geometria. Os primeiros estudos acerca das equações diferenciais se deram nos trabalhos de G. W. Leibniz, I. Newton e dos irmãos Jacob e Johann Bernoulli, os quais determinaram soluções de equações diferenciais de primeira e segunda ordens e aplicaram os seus métodos na descrição de problemas físicos simples ^[1]. No século XVIII, o conhecimento matemático acerca das equações diferenciais deu um salto significativo, devido ao trabalho de renomados matemáticos, dentre os quais se destacam os seguintes nomes: Clariaut, D'Alembert, Lagrange, Laplace, Euler, Gauss, Fourier e Bessel ^[2]. Por exemplo, Euler se destacou por elaborar a técnica do fator integrante para se resolver equações diferenciais de primeira ordem, enquanto Laplace e D'Alembert estudaram técnicas para se resolver equações diferenciais parciais de grande interesse físico, as quais mais tarde receberiam os seus nomes ^[3]. Já no século XIX, Cauchy, Gauss e Liouville também realizaram grandes contribuições no universo das equações diferenciais. Cauchy analisou a convergência de séries infinitas, utilizando-as para encontrar soluções de equações diferenciais com coeficientes variáveis. Gauss evidenciou a importância das funções de variáveis complexas na compreensão das soluções das equações diferenciais aplicadas. Outras figuras importantes no estudo das equações diferenciais foram Picard, Lindelof e Peano, sendo este último o responsável por provar as condições de existência de soluções para equações diferenciais de primeira ordem. Já no século XX, o principal enfoque no campo das equações diferenciais foi quanto aos métodos numéricos, nos quais se destacaram os trabalhos de Carl Runge e Martin Kutta. Desde então, as equações diferenciais passaram a ser aplicadas nas mais diversas áreas do conhecimento, dentre as quais destacam-se a Física, Biologia, Economia, Engenharia e Geologia. Devido ao amplo universo de aplicações, as equações diferenciais tornaram-se umas das principais áreas de estudo da matemática [¹ -³].

Em geral, uma equação diferencial traz uma identidade que envolve derivadas de uma ou mais variáveis dependentes em relação a uma ou mais variáveis independentes. Dessa forma, as equações diferenciais podem ser classificadas em dois tipos, quais sejam: as equações diferenciais ordinárias (EDO's) e as equações diferenciais parciais (EDP's). As equações diferenciais ordinárias envolvem derivadas de uma ou mais variáveis dependentes em relação a uma única variável independente. Já as equações diferenciais parciais envolvem derivadas parciais de uma ou mais variáveis dependentes em relação a mais de uma variável independente. As equações diferenciais são também classificadas por sua ordem e linearidade. A ordem de uma equação diferencial é dada pela ordem da maior derivada da variável dependente que surge na equação ^[3]. Por outro lado, dada uma equação diferencial parcial

(1) $$F\left(x_1,\cdots ,x_n,u,\frac{\partial u}{\partial x_1},\cdots ,\frac{\partial u}{\partial x_n},\frac{{\partial }^{2}u}{\partial x_1^2},\frac{{\partial }^{2}u}{\partial x_1\partial x_2},\cdots \right)=0$$

na variável $u(x_1,\cdots ,x_n)$ dizemos que essa equação é linear se a função $F$ é uma função linear na variável dependente $u$ e nas suas derivadas ^[4]. Encontrar soluções de equações diferenciais constitui um árduo trabalho no qual matemáticos e físicos estão envolvidos. Para esse fim, existem diversos métodos de solução, separados em métodos analíticos e métodos numéricos ^[3].

As equações diferenciais surgem naturalmente na modelagem de diversos tipos de problemas físicos que envolvem sistemas contínuos. Em geral, a modelagem realizada envolve sistemas de equações diferenciais, parciais ou ordinárias, de segunda ordem, muitas vezes relacionados a segunda lei de Newton. Resolver analiticamente as equações obtidas constitui um desafio para os físicos. Nesse arcabouço, a equação de Laplace, a equação de onda, a equação de Schrödinger, a equação de Helmholtz; e a equação de difusão do calor são os exemplos mais corriqueiros de equações diferencias parciais que aparecem na Física. Assim, para solucioná-las, diversos métodos podem ser aplicados a depender das condições iniciais e de contorno estabelecidas [³, ⁴]. Um método semelhante ao que será apresentado aqui é baseado nas chamadas funções de Green ^[5]. Sendo assim, neste trabalho apresentaremos um método alternativo para se resolver analiticamente uma classe de problemas de valor inicial e de contorno para equações diferenciais parciais lineares cujas condições de contorno são essencialmente não--homogêneas.

A grande vantagem do método apresentado neste trabalho está no procedimento de solução que é puramente construtivo. Em geral, problemas de valores de contorno não--homogêneos demandam técnicas de solução relativamente mais complicadas quando comparadas ao simples método da separação de variáveis usualmente utilizado em problemas com equações lineares homogêneas e condições de contorno homogêneas. Logo, neste trabalho iremos mostrar como é possível construir a sua própria transformada integral para determinar a solução do problema desejado. Nesse sentido, o objetivo deste artigo é apresentar, de forma pedagógica, o método das transformadas finitas ^[6] para se resolver problemas de valor inicial e de contorno não-homogêneos para equações diferenciais parciais lineares, de forma geral.

Historicamente, o método das transformadas finitas foi introduzido através das Transformadas Finitas de Fourier [⁷, ⁸] que eram utilizadas para resolver problemas de vibrações em coordenadas Cartesianas. Generalizações desse método foram então levadas para outros tipos de problemas em sistemas de coordenadas distintos, de tal forma que problemas de valor inicial e de contorno envolvendo simetria cilíndrica foram resolvidos utilizando a Transformada Finita de Hankel [⁹], ¹⁰. Além disso, vários problemas de controle envolvendo equações diferenciais lineares, tanto ordinárias quanto parciais, podem também ser resolvidos aplicando a Transformada Finita de Laplace [¹¹, ¹²]. Desta forma, neste texto iremos mostrar como se pode definir uma transformada finita conveniente quando o problema se apresenta em uma forma semelhante a um problema de Sturm--Liouville.

Em particular, por questões didáticas, o problema que será tratado no texto será a equação do calor não--homogênea, devido a sua ampla aplicação na Física e na Engenharia. Lembrando que a técnica aqui apresentada poderá ser aplicada em diversos outros contextos e equações diferenciais distintas contanto que o operador diferencial espacial esteja na forma de um problema de Sturm--Liouville. Observa-se que a técnica da transformada finita aqui apresentada é basicamente uma ótica distinta sobre o método de solução conhecido como expansão por autofunções ^[13]. O método das transformadas finitas busca uma forma de solucionar os problemas de valor inicial e de contorno semelhante ao método das transformadas integrais com espectro contínuo, isto é, simplificando um determinando problema através de relações e propriedades operacionais que a transformada utilizada possui. As Transformadas de Laplace e Fourier ^[13], cujas aplicações são vastamente difundidas na solução de equação diferenciais, são ótimos exemplos das regras operacionais que as transformadas integrais possuem. Essa linguagem diferente pode facilitar bastante a solução do problema, principalmente quando se conhece as propriedades comuns da transformada e também sobre como ela age sobre as derivadas da equação diferencial ^[14], que é um dos objetivos deste trabalho.

Com isso, a apresentação deste trabalho está organizada da seguinte forma. Na seção II será apresentada a equação do calor escalar não-homogênea em um sistema de coordenadas qualquer. Na seção III, trazemos uma breve revisão sobre a teoria de Sturm--Liouville para motivar a construção da transformada. Na seção IV é mostrado como a transformada finita deve ser construída e suas principais propriedades. Na seção V, mostramos utilizando os resultados obtidos na seção IV como é a solução do problema de difusão do calor escalar. Além disso, dois exemplos associados a projetos de pesquisa em andamento na Universidade de Brasília são apresentados e resolvidos usando a técnica apresentada, na seção VI. O primeiro exemplo foi escolhido principalmente pela sua simplicidade, o que facilita o entendimento dos passos relacionados a solução. Nele é possível perceber como o cuidado na hora de se buscar os autovalores do problema pode ser determinante para a solução final. Já o segundo exemplo foi escolhido por ser um sistema acoplado de equações na variável espacial. Neste problema, em particular, a utilização da técnica de desacoplamento do sistema de equações não é possível, o que leva à necessidade de uma forma alternativa de solução. Finalmente, na seção VII são estabelecidas as considerações finais e perspectivas.

2. Equação do calor

Uma das equações diferenciais parciais mais conhecida na Física é a equação do calor, a qual descreve o fluxo de calor em um corpo sólido. A equação do calor é aplicada em muitos contextos, variando desde a construção civil até o estudo do calor dissipado pelo atrito de voos espaciais na reentrada na atmosfera terrestre. Nesta seção, como motivação para a construção da transformada, o nosso enfoque será a apresentação de um problema de difusão de calor unidimensional não-homogêneo considerando condições de contorno também não--homogêneas ^[3].

Nesse sentido, considere, por exemplo, o seguinte problema de valor inicial e contorno relacionado com a difusão de calor em um objeto considerando um sistema de coordenadas espaciais qualquer

(2) $$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}=\kappa {\Delta }_{\circ }u+w_e(r,t),\hfill & (r,t)\in (a,b)\times (0,\infty )\hfill \\ u(r,0)=u_0\left(r\right),\hfill & a\le r\le b\hfill \\ {\alpha }_{1}u(a,t)+{\alpha }_2\frac{\partial u}{\partial r}(a,t)=q_a\left(t\right),\hfill & t\ge 0\hfill \\ {\beta }_{1}u(b,t)+{\beta }_2\frac{\partial u}{\partial r}(b,t)=q_b\left(t\right),\hfill & t\ge 0\hfill \end{array}\right.$$

em que $u(r,t)$ é a temperatura do objeto, $\kappa$ é o coeficiente de difusividade térmica e $w_e(r,t)$ é uma função seccionalmente contínua que representa alguma fonte externa de calor agindo dentro do domínio. Além disso, nas condições de contorno, ${\alpha }_i,$ ${\beta }_i$ são constantes para $i=1,2$ e as funções $q_j\left(t\right)$, com $j=a,b$ são funções seccionalmente contínuas que representam o comportamento da temperatura nos contornos. Além disso, ${\Delta }_{\circ }$ é um operador diferencial relacionado a variável espacial dado por

(3) $${\Delta }_{\circ }u=\frac{1}{W\left(r\right)}\left\{\frac{\partial }{\partial r}\left[P\left(r\right)\frac{\partial u}{\partial r}\right]+Q\left(r\right)u\right\}.$$

Perceba que esse tipo de operador diferencial se assemelha ao operador Laplaciano em coordenadas esféricas ou cilíndricas, por exemplo ^[13]. Neste trabalho vamos considerar que as funções $P\left(r\right)$ e $W\left(r\right)$ são positivas para todo $r\in [a,b]$ e também que $P\left(r\right),P^{\prime }\left(r\right),Q\left(r\right)$ e $W\left(r\right)$ são funções contínuas em todo $r\in [a,b]$

Usualmente, este tipo de problema de valor inicial e contorno é resolvido usando o método da separação de variáveis, isto é, admite-se que a solução da equação diferencial tem a forma $u(r,t)=z\left(t\right)·y\left(r\right)$. No entanto, a simples utilização deste método poderá não ser bem sucedida quando pelo menos uma das funções que definem as condições de contorno $q_a$ e $q_b$ for não-nula ou quando a fonte de calor externa $w_e(r,t)$ for não-nula. Outro método que é aplicável a este caso é o formalismo da função de Green, o qual possui muita visibilidade em livros-texto de física-matemática ^[5]. Nosso enfoque é apresentar uma metodologia diferente e mais intuitiva que possui divulgação menos abrangente na literatura. Desta forma, propomos nesse texto uma metodologia para construir a solução analítica de problemas não--homogêneos e com condições de contorno não--homogêneas como o dado pelo problema (2).

3. Teoria de Sturm--Liouville

Nesta seção, apresentaremos uma breve revisão sobre o problema de Sturm--Liouville. Este formalismo é relevante na solução de muitos problemas físicos, encontrando muitas aplicações nas funções especiais da física-matemática ^[15].

Para este resumo considere que o problema de valor inicial e de contorno (2) não possui fonte de calor externa, isto é, $w_e(r,t)=0$. Ao aplicarmos o método da separação de variáveis nesta versão simplificada do problema chegaremos na seguinte equação diferencial ordinária de segunda ordem

(4) $$\frac{d}{dr}\left[P\left(r\right)\frac{dy}{dr}\right]+Q\left(r\right)y\left(r\right)=\tau W\left(r\right)y\left(r\right)$$

em que $\tau$ é o que chamaremos de um autovalor da equação diferencial. Pode-se provar sem grandes dificuldades que os autovalores relacionados a essa equação são não-positivos, isto é, $\tau =-{\lambda }^2\le 0$, podendo ser estritamente negativo dependendo das condições de contorno [¹⁶, ¹⁷]. Quando as condições de contorno são homogêneas e dadas por

(5) $$\begin{array}{lll}\hfill B_a\left[y\right]:={\alpha }_{1}y\left(a\right)+{\alpha }_2{y}^{\prime }\left(a\right)& \hfill =& \hfill 0\\ \hfill B_b\left[y\right]:={\beta }_{2}y\left(b\right)+{\beta }_2{y}^{\prime }\left(b\right)& \hfill =& \hfill 0\end{array}$$

de tal forma que as constantes ${\alpha }_i$ e ${\beta }_i$, $i=1,2$, satisfazem ${\alpha }_1^2+{\alpha }_2^2>0$ e ${\beta }_1^2+{\beta }_2^2>0$, dizemos que a reunião da equação diferencial dada pela equação (4) com as condições de contorno dadas na equação (5) é um problema de Sturm--Liouville regular, se as funções $P\left(r\right)$ e $W\left(r\right)$ forem positivas para todo $r\in [a,b]$ e que $P\left(r\right),P^{\prime }\left(r\right),Q\left(r\right)$ e $W\left(r\right)$ sejam funções contínuas em todo $r\in [a,b]$ ^[17].

A equação diferencial (4) tem solução única geral dada por

(6) $$y_{\lambda }\left(r\right)=c_{1\lambda }{f}_{\lambda }\left(r\right)+c_{2\lambda }{g}_{\lambda }\left(r\right)$$

em que $f_{\lambda }\left(r\right)$ e $g_{\lambda }\left(r\right)$ são um par de soluções linearmente independentes da equação (4) que estão na base do espaço solução para um autovalor $\lambda$ ^[17]. Os autovalores deste problema são calculados utilizando a solução geral e as condições de contorno. De forma geral, é possível mostrar que os autovalores associados a um problema de Sturm--Liouville formam uma sequência discreta monótona crescente ilimitada [¹⁷, ¹⁸]

(7) $$0\le {\lambda }_1<{\lambda }_2<{\lambda }_3<\cdots <{\lambda }_m<\cdots .$$

Para determinarmos os autovalores considere a autofunção associada ao autovalor ${\lambda }_m$ dada na forma

(8) $$y_m\left(r\right)=c_{1m}{f}_m\left(r\right)+c_{2m}{g}_m\left(r\right),$$

em que $f_m\left(r\right)=f_{{\lambda }_m}\left(r\right)$ e $g_m\left(r\right)=g_{{\lambda }_m}\left(r\right)$. Das condições de contorno temos para $y_m\left(r\right)$ que

(9) $$\begin{array}{l}{\alpha }_{1}ym\left(a\right)+{\alpha }_2{y^{\prime }}_m\left(a\right)=0\\ {\beta }_{1}ym\left(b\right)+{\beta }_2{y^{\prime }}_m\left(b\right)=0\end{array}$$

Considerando $y_m\left(r\right)$ dada pela equação (8), então o sistema ficará

(10) $$\begin{array}{l}{c}_{1m}\left[{\alpha }_1{f}_m\left(a\right)+{\alpha }_2{f^{\prime }}_m\left(a\right)\right]+c_{2m}\left[{\alpha }_{1}gm\left(a\right)+{\alpha }_2{g^{\prime }}_m\left(a\right)\right]=0\\ c_{1m}\left[{\beta }_1{f}_m\left(b\right)+{\beta }_2{f^{\prime }}_m\left(b\right)\right]+c_{2m}\left[{\beta }_{1}gm\left(b\right)+{\beta }_2{g^{\prime }}_m\left(b\right)\right]=0,\end{array}$$

que reescrendo assumirá a forma

(11) $$\begin{array}{l}{c}_{1m}{B}_a\left[f_m\right]+c_{2m}{B}_a\left[g_m\right]=0\\ c_{1m}{B}_b\left[f_m\right]+c_{2m}{B}_b\left[g_m\right]=0.\end{array}$$

Para que esse sistema, nas variáveis $c_{1m}$ e $c_{2m}$, tenha solução não-nula é necessário que o determinante da matriz dos coeficientes seja nulo. A solução não-nula é desejável, caso contrário todas as autofunções do problema serão nulas, o que não é o nosso objetivo encontrar. Desta forma, temos que a equação para calcular os autovalores é dada por

(12) $$B_a\left[f_m\right]B_b\left[g_m\right]-B_a\left[g_m\right]B_b\left[g_m\right]=0.$$

Além disso, para cada autovalor ${\lambda }_m$ estará associada a seguinte autofunção na forma

(13) $$y_m\left(r\right)=B_a\left[f_m\right]g_m\left(r\right)-B_a\left[g_m\right]f_m\left(r\right)$$

que é solução da equação (4) satisfazendo as condições de contorno homogêneas dadas.

Além disso, também é possível mostrar que as autofunções associadas a cada um desses autovalores são duas a duas ortogonais ^[17], isto é, definido o produto interno na forma

(14) $$〈u,v〉={\int }_a^{b}W\left(r\right)u\left(r\right)v\left(r\right)dr,$$

então se $y_n\left(r\right)$ e $y_m\left(r\right)$ são autofunções associadas a autovalores distintos ${\lambda }_n\ne {\lambda }_m$, então o produto interno entre elas satisfaz $〈y_n,y_m〉=0.$ Esse fato é uma consequência do operador ${\Delta }_{\circ }$ ser auto-adjunto [¹⁷, ¹⁸].

O conjunto de autofunções $\left\{y_m\left(r\right):m\in \mathbb{N}\right\}$ na verdade é um conjunto completo no espaço de funções de quadrado integrável ${\mathcal{L}}^2(a,b)=\left\{f:[a,b]\to \mathbb{C}:{\int }_a^b{\left|f\left(x\right)\right|}^{2}dx<\infty \right\}$, isto significa que, o conjunto de autofunções na verdade é uma base de Schauder do espaço de dimensão infinita ${\mathcal{L}}^2(a,b)$ ^[18]. Assim, dada uma função $h\in {\mathcal{L}}^2(a,b)$, então ela pode ser escrita como sendo uma combinação linear das autofunções $y_m$'s na forma

(15) $$h\left(r\right)=\sum _{m=1}^{\infty }{a}_m{y}_m\left(r\right)$$

observando que a convergência é para todo $r\in (a,b)$ [¹⁷, ¹⁸]. É importante observar que essa convergência é uniforme e portanto podemos utilizar a ortogonalidade das autofunções para encontrar os coeficientes $a_m$ ^[17]. Neste caso, podemos facilmente obter que a expansão da função $h\left(r\right)$ no espaço gerado pelo conjunto de autofunções $\left\{y_m\left(r\right):m\in \mathbb{N}\right\}$ é dada na forma

(16) $$h\left(r\right)=\sum _{m=1}^{\infty }\frac{〈f,y_m〉}{〈y_m,y_m〉}{y}_m\left(r\right).$$

Com isso, finalizamos a apresentação da teoria de Sturm-Liouville. Na próxima seção iniciaremos a construção da transformada integral que será utilizada na solução do problema proposto.

4. Construindo a transformada

As transformadas integrais proporcionam poderosas ferramentas para resolver problemas de valores de contorno para equações diferenciais lineares, como o nosso caso. Historicamente, o desenvolvimento e técnicas associados as transformadas integrais estão diretamente relacionados as áreas da matemática aplicada, física e engenharia. Uma transformada integral de uma função $h\left(r\right)$ definida em um intervalo $[a,b]$ é o operador linear definido por

(17) $$\mathcal{T}\left[h\left(r\right)\right]={\int }_a^{b}K(r,\lambda )h\left(r\right)dr$$

em que $K(r,\lambda )$ é chamado de núcleo da transformada, ou kernel ^[14]. Desejamos construir uma transformada que possa ser aplicada, por exemplo, na equação diferencial homogênea

(18) $$\frac{\partial u}{\partial t}=\kappa {\Delta }_{\circ }u$$

e que ofereça uma versão transformada que seja conveniente para que possamos encontrar de forma mais simples a solução do problema. Como sabemos, pela discussão das seções anteriores, que o operador diferencial dado pela equação (3) está de alguma forma associado com o problema de Sturm--Liouville descrito na seção 3, então usaremos as soluções do problema de Sturm--Liouville homogêneo como base para a nossa transformada. Desta forma, definimos a transformada ${\mathcal{I}}_m[h]$ como sendo

(19) $${\mathcal{T}}_m\left[h\right]=\tilde{h}\left({\lambda }_m\right)={\int }_a^{b}W\left(r\right)h\left(r\right)y_m\left(r\right)dr,$$

em que $y_m\left(r\right)=B_a\left[f_m\right]g_m\left(r\right)-B_a\left[g_m\right]f_m\left(r\right)$ e os autovalores ${\lambda }_m$ são calculados usando a relação dada pela equação (12). Neste caso, o núcleo da nossa transformada é dado por $K\left(r,{\lambda }_m\right)=W\left(r\right)y_m\left(r\right)$, em que a função peso $W\left(r\right)$ é dada pela equação (3).

Por definição essa transformada é linear, isto é, dados $\eta ,\mu \in \mathbb{R}$ constantes, tem-se ${\mathcal{T}}_m\left[\xi u\left(r\right)+\zeta v\left(r\right)\right]=\xi {\mathcal{T}}_m\left[u\left(r\right)\right]+\zeta {\mathcal{T}}_m\left[v\left(r\right)\right].$ Com o objetivo de aplicar essa transformada na equação diferencial é necessário saber como ela se comporta no operador diferencial ${\Delta }_{\circ }h.$ Assim, integrando por partes algumas vezes, temos que a transformada no operador diferencial satisfará

(20) $$\begin{array}{ccc}& & {\mathcal{T}}_m\left[{\Delta }_{\circ }h\right]={\int }_a^{b}W\left(r\right)\left({\Delta }_{\circ }h\right)y_m\left(r\right)dr\hfill \\ \\ & & ={\int }_a^b\left\{\frac{d}{dr}\left[P\left(r\right)\frac{dh}{dr}\right]+Q\left(r\right)h\left(r\right)\right\}{y}_m\left(r\right)dr\hfill \\ & & ={\left[P\left(r\right)h^{\prime }\left(r\right)y_m\left(r\right)\right]}_a^b-{\int }_a^{b}P\left(r\right)\frac{dh}{dr}\frac{d{y}_m}{dr}dr\hfill \\ & & +{\int }_a^{b}Q\left(r\right)h\left(r\right)y_m\left(r\right)dr\hfill \\ & & ={\left[P\left(r\right)h^{\prime }\left(r\right)y_m\left(r\right)\right]}_a^b-{\left[P\left(r\right)h\left(r\right)y_m^{\prime }\left(r\right)\right]}_a^b\hfill \\ & & +{\int }_a^{b}h\left(r\right)\left\{\frac{d}{dr}\left[P\left(r\right)\frac{d{y}_m}{dr}\right]+Q\left(r\right)y_m\left(r\right)\right\}dr\hfill \\ & & ={\left[P\left(r\right)h^{\prime }\left(r\right)y_m\left(r\right)\right]}_a^b-{\left[P\left(r\right)h\left(r\right)y_m^{\prime }\left(r\right)\right]}_a^b\hfill \\ & & +{\int }_a^{b}h\left(r\right)\left(-{\lambda }_m^{2}W\left(r\right)y_m\left(r\right)\right)dr\hfill \\ & & =-{\lambda }_m^2{\mathcal{T}}_m\left[h\right]+{\left[P{h}^{\prime }{y}_m\right]}_a^b-{\left[Ph{y}_m^{\prime }\right]}_a^b.\hfill \end{array}$$

Considerando as condições de contorno, podemos ainda simplificar um pouco mais a relação acima para obter

(21) $${\mathcal{T}}_m\left[{\Delta }_{\circ }h\right]=-{\lambda }_m^2{\mathcal{T}}_m\left[h\right]+\frac{1}{{\beta }_2}{\left[P\left(r\right)y_m\left(r\right)B_r\left[h\right]\right]}_a^b,$$

considerando que ${\beta }_2\ne 0$.

Perceba que a definição da transformada discreta nada mais é que que o produto interno entre $y_m$ e $h$, isto é ${\mathcal{T}}_m\left[h\right]=〈h,y_m〉$. Desta forma, a transformada inversa é dada pela expansão

(22) $$h\left(r\right)={\mathcal{T}}_m^{-1}\left[\tilde{h}\left({\lambda }_m\right)\right]=\sum _{m=1}^{\infty }\frac{{\mathcal{T}}_m\left[h\right]}{{\mathcal{T}}_m\left[y_m\right]}{y}_m\left(r\right).$$

Essa última equação será utilizada a partir da próxima seção para solucionarmos o problema da condução do calor. Perceba que a transformada finita nada mais é que uma visão diferente sobre a expansão por autofunções. Veremos nos exemplos como esse conceito se faz conveniente na hora de solucionar os problemas de valor inicial e contorno.

5. Resolvendo o problema

Podemos agora resolver o nosso problema aplicando a transformada definida na Seção 4. Considerando que $\tilde{u}({\lambda }_m,t)={\mathcal{T}}_m\left[u(r,t)\right]$ e ${\tilde{w}}_e({\lambda }_m,t)={\mathcal{T}}_m\left[w_e(r,t)\right]$, então aplicando a transformada na equação diferencial parcial obtemos a seguinte equação diferencial ordinária de primeira ordem na variável $t$

(23) $$\frac{\partial \tilde{u}}{\partial t}=-\kappa {\lambda }_m^2\tilde{u}+\frac{1}{{\beta }_2}{\left[P\left(r\right)y_m\left(r\right)q_r\left(t\right)\right]}_a^b+{\tilde{w}}_e({\lambda }_m,t).$$

Essa equação pode ser facilmente resolvida usando o método do fator integrante. Para tal considerando $\mu \left(t\right)=e^{\kappa {\lambda }_m^{2}t}$ e integrando a equação obtemos

(24) $$\begin{array}{ccc}\hfill \tilde{u}({\lambda }_m,t)& =& \tilde{c}\left({\lambda }_m\right)e^{-\kappa {\lambda }_m^{2}t}+\frac{{\left[P\left(r\right)y_m\left(r\right)\right]}_a^b}{{\beta }_2}\times \hfill \\ & & \times {\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\left[q_r\left(s\right)\right]}_a^{b}ds\hfill \\ & & +{\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\tilde{w}}_e({\lambda }_m,s)ds,\hfill \end{array}$$

em que $\tilde{c}\left({\lambda }_m\right)$ é uma constante de integração. Se considerarmos que $\tilde{u}({\lambda }_m,0)$ é a transformada da condição inicial $\tilde{u}({\lambda }_m,0)={\tilde{u}}_0\left({\lambda }_m\right)={\mathcal{T}}_m\left[u_0\left(r\right)\right]$, então podemos determinar a constante de integração e a solução da equação será dada por

(25) $$\begin{array}{ccc}\hfill \tilde{u}({\lambda }_m,t)& =& {\tilde{u}}_0\left({\lambda }_m\right)e^{-\kappa {\lambda }_m^{2}t}+\frac{{\left[P\left(r\right)y_m\left(r\right)\right]}_a^b}{{\beta }_2}\times \hfill \\ & & \times {\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\left[q_r\left(s\right)\right]}_a^{b}ds\hfill \\ & & +{\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\tilde{w}}_e({\lambda }_m,s)ds.\hfill \end{array}$$

Finalmente, podemos escrever a solução final do problema de valor inicial e de contorno que é dada por

(26) $$\begin{array}{ccc}\hfill u(r,t)& =& \sum _{m=1}^{\infty }\left\{{\tilde{u}}_0\left({\lambda }_m\right)e^{-\kappa {\lambda }_m^{2}t}+\frac{{\left[P\left(r\right)y_m\left(r\right)\right]}_a^b}{{\beta }_2}\times \right.\hfill \\ & & \times {\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\left[q_r\left(s\right)\right]}_a^{b}ds\hfill \\ & & +\left.{\int }_0^t{e}^{-\kappa {\lambda }_m^2(t-s)}{\tilde{w}}_e({\lambda }_m,s)ds\right\}{Y}_m\left(r\right)\hfill \end{array}$$

em que $Y_m\left(r\right)=y_m\left(r\right)/〈y_m,y_m〉$ é a autofunção normalizada do problema de Sturm-Liouville caracterizado pelas equações (4)-(5).

6. Exemplos

Nesta seção iremos determinar as soluções de dois problemas distintos usando a técnica da transformada finita. Os problemas foram escolhidos por estarem em sistemas de coordenadas distintos e também por possuírem condições de contorno distintas. Desta forma, diferentes sequências de autofunções e autovalores surgirão para cada um dos casos.

O primeiro problema, apresentado na seção 6.1, está associado ao estudo do escudo térmico para um motor foguete que está sendo estudado e construído na Universidade de Brasília ^[19]. Este é um problema que não possui fonte de calor externa na equação diferencial, mas possui condições de contorno não--homogêneas e também simetria cilíndrica. Ele foi escolhido para ser descrito nesse trabalho mais pela sua qualidade didática na descrição da técnica do que pelo seu grau de dificuldade.

Por outro lado, na seção 6.2, trataremos da solução de um sistema de equações relacionado ao problema da ablação hepática por rádiofrequência ^[20], cujo estudo na Universidade de Brasília está associado a construção de um protótipo para o equipamento médico. Este problema é não--homogêneo, possui simetria esférica e é formado por um sistema de equações acopladas ^[21]. Neste caso, a matriz que realiza o acoplamento não é diagonalizável levando a uma necessidade de encontrar um método distinto para a resolução do problema. Portanto, o método da transformada finita descrito nesse trabalho se mostra uma ferramenta interessante para contornar essa situação.

6.1 Equação do calor escalar em coordenadas polares

Para fixar o entendimento do método, neste primeiro exemplo vamos considerar o seguinte problema de valor inicial e de contorno

(6.1.1) $$\left\{\begin{array}{cc}\frac{1}{\alpha }\frac{\partial T}{\partial t}={\Delta }_{c}T,\hfill & (r,t)\in (a,b)\times (0,\infty )\hfill \\ T(r,0)=f_0\left(r\right),\hfill & 0<a\le r\le b\hfill \\ \frac{\partial T}{\partial r}(a,t)=q_a,\hfill & t\ge 0\hfill \\ \frac{\partial T}{\partial r}(b,t)=q_b,\hfill & t\ge 0\hfill \end{array}\right..$$

Este é um problema de difusão de calor em uma geometria cilíndrica na forma de um anel, como pode ser visto na figura (1). Conforme citado previamente, esse problema está relacionado com o estudo do isolamento térmico em um motor foguete, isto é deseja-se compreender o comportamento da temperatura na região entre $0<a\le r\le b$ sujeito ao fluxo de calor intenso gerado pelo motor foguete que passa pela região entre $0\le r<a$ ^[19]. Neste caso, $\alpha$ é o coeficiente de difusividade térmica do material isolante e $T(r,t)$ a temperatura que o material isolante terá em função de $(r,t)$.

Figura 1
Domínio em que se deseja encontrar a distribuição de temperatura em função do tempo. O material do anel é um isolante térmico enquanto a região compreendida entre $0\le r\le a$ é uma região de fluxo de calor intenso gerado pelo motor foguete ^[19]

Aqui a condição inicial é dada por uma temperatura $T(r,0)=f_0\left(r\right)$ e as condições de contorno são dadas por ${\displaystyle {\left.\frac{\partial T}{\partial r}\right|}_a=q_a}$ e ${\displaystyle {\left.\frac{\partial T}{\partial r}\right|}_b=q_b}$, em que $q_a$ e $q_b$ são os fluxos de calor constantes nas paredes interna e externa do cilindro.

Com o objetivo de facilitar a compreensão do método da transformada finita vamos dividir a construção da solução do problema (6.1.1) em algumas etapas. Começaremos estudando os autovalores do problema de Sturm-Liouville associado ao problema (6.1.1) desejado. Em seguida, as autofunções deste mesmo problema de Sturm-Liouville serão determinadas, pois elas são fundamentais para a construção da transformada. Na próxima etapa a transformada discreta será definida, para finalmente ser aplicada ao problema (6.1.1) e determinarmos a expressão final da solução.

6.1.1 Tipos de autovalores

O primeiro passo para obter a solução geral do problema (6.1.1) é entender o caráter dos autovalores do problema de Sturm--Liovuille associado a ele. Esta etapa é importante, pois a transformada discreta é construída baseada nas autofunções associadas a esses autovalores. Desta forma, uma análise incorreta do conjunto de autovalores pode levar a uma solução equivocada ao final do processo.

Neste caso, temos que o operador diferencial ${\Delta }_c$ é o operador Laplaciano em simetria cilíndrica que é dado por

(6.1.2) $${\Delta }_{c}h=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial h}{\partial r}\right),$$

isto é, as funções $P\left(r\right),Q\left(r\right)$ e $W\left(r\right)$ são dadas por $P\left(r\right)=r$, $Q\left(r\right)=0$ e $W\left(r\right)=r$. Sendo assim, a solução do problema (6.1.1) está ligada, conforme visto nas Seções 3 e 4, à solução da equação de Bessel de ordem zero abaixo

(6.1.3) $${\left[r{R}^{\prime }\left(r\right)\right]}^{\prime }=\tau rR\left(r\right)$$

sujeito as condições de contorno homogêneas $R^{\prime }\left(a\right)=R^{\prime }\left(b\right)=0$. Aqui é importante ressaltar que o domínio do problema satisfaz $0<a<b$ e portanto, nesse intervalo, o problema de Sturm--Liouville é regular.

Desta forma, multiplicando a equação diferencial por $R\left(r\right)$ e integrando, tem-se

(6.1.4) $${\int }_a^{b}R\left(r\right){\left[r{R}^{\prime }\left(r\right)\right]}^{\prime }dr=\tau {\int }_a^{b}rR{\left(r\right)}^{2}dr.$$

Agora, integrando por partes o lado esquerdo da equação resulta em

(6.1.5) $${\left[rR\left(r\right)R^{\prime }\left(r\right)\right]}_a^b-{\int }_a^{b}r{R}^{\prime }{\left(r\right)}^{2}dr=\tau {\int }_a^{b}rR{\left(r\right)}^{2}dr.$$

Como as condições de contorno são $R^{\prime }\left(a\right)=R^{\prime }\left(b\right)=0$, tem-se que ${\left[rR\left(r\right)R^{\prime }\left(r\right)\right]}_a^b=0.$ Portanto, os autovalores do problema de Sturm--Liouville associado satisfazem

(6.1.6) $$\tau =-\frac{{\int }_a^{b}r{R}^{\prime }{\left(r\right)}^{2}dr}{{\int }_a^{b}rR{\left(r\right)}^{2}dr}\le 0$$

e conclui-se que os autovalores devem ser não-positivos para esse tipo de problema, conforme já havíamos citado anteriormente.

6.1.2 Autovalores, autofunções e ortogonalidade

Sabendo que os autovalores são não-positivos, será analisado inicialmente como são calculados os autovalores estritamente negativos $\lambda <0$, que serão escritos a partir de agora como $\tau =-{\lambda }^2$. Neste caso, a equação (6.1.3) pode ser escrita como

(6.1.7) $${\left[r{R}^{\prime }\left(r\right)\right]}^{\prime }=-{\lambda }^{2}rR\left(r\right).$$

Conforme observado anteriormente, essa é uma equação de Bessel de ordem zero que possui solução geral dada por

(6.1.8) $$R\left(r\right)=c_1{J}_0\left(\lambda r\right)+c_2{Y}_0\left(\lambda r\right),$$

em que $J_0\left(\lambda r\right)$ e $Y_0\left(\lambda r\right)$ são as funções de Bessel de primeira e segunda espécie de ordem zero ^[16]. Para determinar os autovalores $\lambda$ deste problema utiliza-se as duas condições de contorno homogêneas $R^{\prime }\left(a\right)=R^{\prime }\left(b\right)=0$ que resulta no seguinte sistema linear

(6.1.9) $$\begin{array}{ccc}\hfill c_1{J}_1\left(\lambda a\right)+c_2{Y}_1\left(\lambda a\right)& =& 0\hfill \\ \hfill c_1{J}_1\left(\lambda b\right)+c_2{Y}_1\left(\lambda b\right)& =& 0,\hfill \end{array}$$

lembrando que as derivadas das funções de Bessel de ordem zero satisfazem ${\displaystyle {\left[J_0\left(\lambda r\right)\right]}^{\prime }=-\lambda J_1\left(\lambda r\right)}$ e ${\displaystyle {\left[Y_0\left(\lambda r\right)\right]}^{\prime }=-\lambda Y_1\left(\lambda r\right)}$ ^[16] que são as funções de Bessel de primeira e segunda espécie de ordem 1. Da segunda equação em (6.1.9), fazendo

(6.1.10) $$c_1=-\frac{c_2{Y}_1\left(\lambda b\right)}{J_1\left(\lambda b\right)},$$

tem-se que os autovalores estritamente negativos são calculados através da relação

(6.1.11) $$J_1\left(\lambda a\right)Y_1\left(\lambda b\right)-J_1\left(\lambda b\right)Y_1\left(\lambda a\right)=0,$$

conforme esperado pela equação (12). Importante observar que os autovalores nesse caso só podem ser calculados numericamente. Para tal pode ser usado que os autovalores de um problema de Sturm--Liouville são assintoticamente dados por ^[22]

(6.1.12) $${\lambda }_m\sim \frac{m\pi }{b-a}$$

em conjunto com o método de Newton--Raphson ^[23] para determinar as raízes da equação (6.1.11).

Associados a esses autovalores estritamente negativos tem-se também a sequência de autofunções

(6.1.13) $$R_m\left(r\right)=J_0\left({\lambda }_{m}r\right)Y_1\left({\lambda }_{m}b\right)-J_1\left({\lambda }_{m}b\right)Y_0\left({\lambda }_{m}r\right),m\ge 1,$$

que são consequência da equação (13).

Por outro lado, temos que o problema pode admitir o autovalor nulo, isto é ${\lambda }_0=0$. Para este caso a equação diferencial do problema se reduz a

(6.1.14) $${\left[r{R}_0^{\prime }\left(r\right)\right]}^{\prime }=0.$$

Portanto, integrando os dois lados tem-se

(6.1.15) $$r{R}_0^{\prime }\left(r\right)=d_1$$

e, finalmente,

(6.1.16) $$R_0\left(r\right)=d_{1}lnr+d_2.$$

Para que essa função satisfaça as condições de contorno dadas é necessário que $d_1=0$. Desta forma, existe uma autofunção associada ao autovalor nulo é dada por $R_0\left(r\right)=1$ que também é solução da equação diferencial (6.1.3) e satisfaz as condições de contorno $R^{\prime }\left(a\right)=R^{\prime }\left(b\right)=0$.

Para verificarmos a ortogonalidade destas autofunções obtidas definimos o produto interno entre duas funções $u\left(r\right)$ e $v\left(r\right)$ como sendo

(6.1.17) $$〈u,v〉={\int }_a^{b}ru\left(r\right)v\left(r\right)dr,$$

pois $W\left(r\right)=r$. Desta forma, observa-se de imediato que a autofunção $R_0\left(r\right)=1,$ por exemplo, associada ao autovalor ${\lambda }_0=0$ é ortogonal a todas as demais $R_m\left(r\right)$, pois

(6.1.18) $$\begin{array}{ccc}& & {\int }_a^{b}r{R}_0\left(r\right)R_m\left(r\right)dr={\int }_a^{b}r·1·R_m\left(r\right)dr\hfill \\ & & ={\int }_a^{b}r\left(J_0\left({\lambda }_{m}r\right)Y_1\left({\lambda }_{m}b\right)-J_1\left({\lambda }_{m}b\right)Y_0\left({\lambda }_{m}r\right)\right)dr\hfill \\ & & ={\left[\frac{r}{{\lambda }_m}\left(J_1\left({\lambda }_{m}r\right)Y_1\left({\lambda }_{m}b\right)-J_1\left({\lambda }_{m}b\right)Y_1\left({\lambda }_{m}r\right)\right)\right]}_a^b\hfill \\ & & =0,\hfill \end{array}$$

considerando que ${\displaystyle \int r{J}_0\left({\lambda }_{m}r\right)dr=\frac{r}{{\lambda }_m}{J}_1\left({\lambda }_{m}r\right)}+C$ e ${\displaystyle \int r{Y}_0\left({\lambda }_{m}r\right)dr=\frac{r}{{\lambda }_m}{Y}_1\left({\lambda }_{m}r\right)}+C$ ^[16].

Neste caso, a condição inicial $T(r,0)=f_0\left(r\right)$ pode ser expandida através da série de funções ortogonais

(6.1.19) $$f_0\left(r\right)=c_0{R}_0\left(r\right)+c_1{R}_1\left(r\right)+\cdots +c_n{R}_n\left(r\right)+\cdots .$$

Como as funções são duas a duas ortogonais, pode-se encontrar os coeficientes da expansão que são dados por

(6.1.20) $$〈R_0,f_0〉=c_0〈R_0,R_0〉$$

e

(6.1.21) $$〈R_m,f_0〉=c_m〈R_m,R_m〉,m>0.$$

No primeiro caso,

(6.1.22) $$\begin{array}{ccc}\hfill c_0& =& \frac{〈R_m,f_0〉}{〈R_0,R_0〉}=\frac{{\int }_a^{b}r{f}_0\left(r\right)·1·dr}{{\int }_a^{b}r·{\left(1\right)}^{2}dr}\hfill \\ & =& \frac{2}{b^2-a^2}{\int }_a^{b}r{f}_0\left(r\right)dr\hfill \end{array}$$

e nos demais, $m\ge 1$,

(6.1.23) $$c_m=\frac{〈R_m,f_0〉}{〈R_m,R_m〉}.$$

Portanto, a expansão da condição inicial é dada na forma

(6.1.24) $$f_0\left(r\right)=\frac{2}{b^2-a^2}{\int }_a^{b}r{f}_0\left(r\right)dr+\sum _{m=1}^{\infty }〈R_m,f_0〉N_m\left(r\right),$$

em que $N_m\left(r\right)=R_m/〈R_m,R_m〉$ é o conjunto das autofunções normalizadas. Importante observar que essa expansão será fundamental para a obtenção da solução final.

6.1.3. Transformada discreta

Agora podemos definir a transformada discreta para o problema (6.1.1) que é dada por

(6.1.25) $${\mathcal{H}}_m\left[h\right]={\int }_a^{b}rh\left(r\right)R_m\left(r\right)dr,$$

em que o núcleo da transformada é dado pela equação (6.1.13) e o autovalores ${\lambda }_m$'s são calculados através da relação dada pela equação (6.1.11). Para o caso $m=0$, temos que $R_0\left(r\right)=1$ e ${\lambda }_0=0$. Assim, temos que a transformada inversa é dada pela expansão

(6.1.26) $$h\left(r\right)=\frac{2}{b^2-a^2}{\int }_a^{b}rh\left(r\right)dr+\sum _{m=1}^{\infty }\frac{{\mathcal{H}}_m\left[h\right]}{〈R_m,R_m〉}{R}_m\left(r\right).$$

Da equação (21), temos que no operador diferencial, para $m>0$, a transformada se comporta na forma

(6.1.27) $${\mathcal{H}}_m\left[{\Delta }_{c}h\right]=-{\lambda }_m^2{\mathcal{H}}_m\left[h\right]+{\left[r\frac{\partial h}{\partial r}{R}_m\right]}_a^b,m>0.$$

No caso do autovalor nulo, é imediato mostrar que para ${\lambda }_0=0,$ a transformada se comporta no operador diferencial da forma

(6.1.28) $${\mathcal{H}}_0\left[{\Delta }_{c}h\right]=b{\left.\frac{\partial h}{\partial r}\right|}_b-a{\left.\frac{\partial h}{\partial r}\right|}_a.$$

6.1.4 Transformando a equação diferencial

Finalmente, podemos aplicar a transformada definida pela equação (6.1.25) na equação do calor dada no problema (6.1.1). Neste caso, tem-se

(6.1.29) $$\frac{1}{\alpha }{\tilde{T}}_t({\lambda }_m,t)=-{\lambda }_m^2\tilde{T}({\lambda }_m,t)+{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b,m\ge 0,$$

considerando que $\tilde{T}({\lambda }_m,t)={\mathcal{H}}_m\left[T(r,t)\right]$. Reescrevendo essa equação fica na forma

(6.1.30) $${\tilde{T}}_t({\lambda }_m,t)+\alpha {\lambda }_m^2\tilde{T}({\lambda }_m,t)=\alpha {\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b.$$

Para $m>0$, temos que a equação acima é uma equação diferencial ordinária de primeira ordem na variável $t$, desta forma multiplicando-a pelo fator integrante $\mu \left(t\right)={\displaystyle e^{\alpha {\lambda }_m^{2}t}}$, tem-se que

(6.1.31) $$\frac{d}{dt}\left(\tilde{T}({\lambda }_m,t)e^{\alpha {\tau }_m^{2}t}\right)=\alpha {\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b{e}^{\alpha {\lambda }_m^{2}t}.$$

Considerando que as condições de contorno $T_r(b,t)$ e $T_r(a,t)$ são constantes e integrando a equação acima tem-se que

(6.1.32) $$\tilde{T}({\lambda }_m,t)e^{\alpha {\lambda }_m^{2}t}=\frac{1}{{\lambda }_m^2}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b{e}^{\alpha {\lambda }_m^{2}t}+\tilde{h}\left({\lambda }_m\right),$$

em que $\tilde{h}\left({\lambda }_m\right)$ é uma constante de integração. Portanto,

(6.1.33) $$\tilde{T}({\lambda }_m,t)=\frac{1}{{\lambda }_m^2}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b+\tilde{h}\left({\lambda }_m\right)e^{-\alpha {\lambda }_m^{2}t},m>0.$$

Agora, por outro lado, para o autovalor nulo, isto é $m=0$, a equação se reduz a

(6.1.34) $${\tilde{T}}_t({\lambda }_0,t)=\alpha {\left[r\frac{\partial T}{\partial r}\right]}_a^b$$

e consequentemente, sua solução é dada por

(6.1.35) $$\tilde{T}({\lambda }_0,t)=\alpha {\left[r\frac{\partial T}{\partial r}\right]}_a^{b}t+\tilde{g}\left({\lambda }_m\right),$$

em que $\tilde{g}\left({\lambda }_m\right)$ novamente é uma constante de integração. Portanto, tem-se que a solução é dada por

(6.1.36) $$\begin{array}{ccc}\hfill & & T(r,t)=\sum _{m=0}^{\infty }\tilde{T}({\lambda }_m,t)N_m\left(r\right)\hfill \\ \hfill & & =\left\{\alpha {\left[r\frac{\partial T}{\partial r}\right]}_a^{b}t+\tilde{g}\left({\lambda }_m\right)\right\}{N}_0\left(r\right)\hfill \\ \hfill & & +\sum _{m=1}^{\infty }\left\{\frac{1}{{\lambda }_m^2}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b+\tilde{h}\left({\lambda }_m\right)e^{-\alpha {\lambda }_m^{2}t}\right\}{N}_m\left(r\right)\hfill \\ \hfill & & =\frac{2}{b^2-a^2}\left\{\alpha {\left[r\frac{\partial T}{\partial r}\right]}_a^{b}t+\tilde{g}\left({\lambda }_m\right)\right\}\hfill \\ \hfill & & +\sum _{m=1}^{\infty }\left\{\frac{1}{{\lambda }_m^2}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b+\tilde{h}\left({\lambda }_m\right)e^{-\alpha {\lambda }_m^{2}t}\right\}\hfill \\ \hfill & & \times N_m\left(r\right).\hfill \end{array}$$

No instante inicial, $t=0$, a expansão da temperatura inicial $f_0\left(r\right)$ é dada por

(6.1.37) $$f_0\left(r\right)=\frac{2}{b^2-a^2}{\int }_a^{b}r{f}_0\left(r\right)dr+\sum _{m=1}^{\infty }〈R_m,f_0〉N_m\left(r\right).$$

Comparando a série da equação (6.1.37) com a série obtida anteriormente para a distribuição de temperatura, equação (6.1.36), pode-se determinar quem são as constantes de integração $\tilde{g}\left({\lambda }_m\right)$ e $\tilde{h}\left({\lambda }_m\right)$, que são dadas por

(6.1.38) $$\begin{array}{ccc}\hfill \tilde{g}\left({\lambda }_m\right)& =& {\int }_a^{b}r{f}_0\left(r\right)dr\hfill \end{array}$$

(6.1.39) $$\begin{array}{ccc}\hfill \tilde{h}\left({\lambda }_m\right)& =& 〈R_m,f_0〉-\frac{1}{{\lambda }_m}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b.\hfill \end{array}$$

Finalmente, podemos escrever a solução final que é dada por

(6.1.40) $$\begin{array}{ccc}\hfill T(r,t)& =& \frac{2}{b^2-a^2}\left\{\alpha {\left[r\frac{\partial T}{\partial r}\right]}_a^{b}t+{\int }_a^{b}r{f}_0\left(r\right)dr\right\}\hfill \\ & +& \sum _{m=1}^{\infty }\frac{1}{{\lambda }_m^2}{\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b\left(1-e^{-\alpha {\lambda }_m^{2}t}\right)N_m\left(r\right)\hfill \\ & +& \sum _{m=1}^{\infty }〈R_m,f_0〉e^{-\alpha {\lambda }_m^{2}t}{N}_m\left(r\right).\hfill \end{array}$$

Além disso, é possível mostrar com algum trabalho usando as relações de produto cruzado e Wronskiano das funções de Bessel que ^[24]

(6.1.41) $$〈R_m,R_m〉=\frac{2}{{\pi }^2{\lambda }_m^2}\left[1-\frac{J_1^2\left({\lambda }_{m}b\right)}{J_1^2\left({\lambda }_{m}a\right)}\right]$$

e

(6.1.42) $${\left[r\frac{\partial T}{\partial r}{R}_m\right]}_a^b=\frac{2}{\pi {\lambda }_m}\left[{\left.\frac{\partial T}{\partial r}\right|}_a\frac{J_1\left({\lambda }_{m}b\right)}{J_1\left({\lambda }_{m}a\right)}-{\left.\frac{\partial T}{\partial r}\right|}_b\right].$$

Portanto, a solução final do problema (6.1.1) pode ser reescrita na forma

(6.1.43) $$\begin{array}{ccc}& & T(r,t)=\frac{2\alpha t}{b^2-a^2}\left[b{q}_b-a{q}_a\right]+\frac{2}{b^2-a^2}{\int }_a^{b}r{f}_0\left(r\right)dr\hfill \\ \hfill \\ & & +\pi \sum _{m=1}^{\infty }\frac{\left[q_a\frac{J_1\left({\lambda }_{m}b\right)}{J_1\left({\lambda }_{m}a\right)}-q_b\right]J_1^2\left({\lambda }_{m}a\right)\left(1-e^{-\alpha {\lambda }_m^{2}t}\right)R_m\left(r\right)}{{\lambda }_m\left[J_1^2\left({\lambda }_{m}a\right)-J_1^2\left({\lambda }_{m}b\right)\right]}\hfill \\ \hfill \\ & & +\frac{{\pi }^2}{2}\sum _{m=1}^{\infty }〈R_m,f_0〉\frac{{\lambda }_m^2{J}_1^2\left({\lambda }_{m}a\right)e^{-\alpha {\lambda }_m^{2}t}{R}_m\left(r\right)}{\left[J_1^2\left({\lambda }_{m}a\right)-J_1^2\left({\lambda }_{m}b\right)\right]},\hfill \end{array}$$

em que ${\lambda }_m$ são os autovalores e $R_m\left(r\right)$ são as autofunções normalizadas calculadas, respectivamente, através das equações (6.1.11) e (6.1.13). Na figura (2), podemos ver a solução do problema (6.1.1), nas extremidades do domínio, para dois diferentes tipos de condições de contorno.

Figura 2
Solução dada pela equação (6.1.43) em duas posições distintas do domínio e para duas condições de contorno distintas: (i) Curvas em azul e vermelho referem-se ao problema (6.1.1) com temperatura inicial constante $f_0\left(r\right)=30{}^{\circ }$ condições de contorno $q_a=-31428.57{}^{\circ }$ e $q_b=0$; (ii) Curvas em preto e verde referem-se ao problema (6.1.1) com temperatura inicial constante $f_0\left(r\right)=30{}^{\circ }$, condições de contorno $q_a=-31428.57[H(t-16)-H(t-6)]{}^{\circ }$ C/m e $q_b=0$, em que $H\left(t\right)$ é a função degrau unitário. Em ambos os casos a condutividade térmica é dada por $\alpha =1.77041286\times {10}^{-7}$ m${}^2$/s, o raio interno $a=38.57$ mm e o raio externo $b=43.57$ mm ^[19].

Ressaltamos que a função $T(r,t)$ encontrada descreve a distribuição de temperatura dentro do anel, mostrado na figura (1), ao longo da distância radial $r$ e dependente do tempo $t$. Além disso, observamos que nesse exemplo a convergência da série dada pela equação (6.1.43) é para todo $r\in [a,b]$ e $t\ge 0$.

6.2. Sistema de Equações Radiais

Vamos tratar nesse segundo exemplo de um modelo matemático utilizado para o estudo dinâmico da propagação de calor na ablação hepática por radiofrequência ^[21]. Neste modelo o tecido hepático humano é considerado como sendo um meio poroso composto por duas fases, tecido e sangue. No interior do fígado uma fonte de calor esférica de raio $0<r<a$ é acionada levando ao aumento de temperatura da região que a envolve. Desta forma, deseja-se estudar como a temperatura do tecido hepático varia com o tempo nas proximidades da fonte de calor no intervalo $[a,b]$, com $a>0$. Neste caso, considera-se a existência de simetria esférica e que a temperatura de cada fase é então modelada pelo seguinte sistema de equações diferenciais

(6.2.1) $$\begin{array}{l}\delta {\rho }_t{c}_t\frac{\partial w_t}{\partial t}=\delta k_t{\Delta }_e{w}_t-h_{tb}\left(w_t-w_b\right)+\delta qt+\delta \frac{v_t}{r^4}\\ ∊{\rho }_b{c}_b\frac{\partial w_b}{\partial t}=∊k_b{\Delta }_e{w}_b+h_{tb}\left(w_t-w_b\right)+∊q_b+∊\frac{v_b}{r^4},\end{array}$$

em que $w_i$ é a temperatura da fase, ${\rho }_i$ é densidade, $c_i$ o calor específico, $k_i$ a difusividade térmica, $h_{tb}$ é o coeficiente de transferência de calor interfacial, $q_i$ é alguma fonte constante de geração interna de calor, ${\nu }_i$ está relacionado a condutividade elétrica da fase, com $i=t$ para o tecido e $i=b$ para o sangue. Além disso, temos que a porosidade é dada pelas constantes $\delta$ e $\varepsilon$ que devem satisfazer $\delta +\varepsilon =1$. É importante aqui observar que as constantes físicas associadas a cada fase serão consideradas diferentes mesmo que em alguns casos elas possam ser próximas ^[25].

O modelo dado pelas equações em (6.2.1) é uma generalização em coordenadas esféricas dos modelos encontrados nas referências [²¹, ²⁶] que considera apenas o problema em regime permanente sem os efeitos elétricos associados. A solução será baseada no problema de valor inicial e de contorno que considera

(6.2.2) $$\left\{\begin{array}{cc}{w}_t(a,t)=w_b(a,t)=w_0\left(t\right),\hfill & t\ge 0\hfill \\ w_t(b,t)=w_b(b,t)=w_c,\hfill & t\ge 0\hfill \\ w_t(r,0)=w_b(r,0)=f\left(r\right),\hfill & r\in [a,b],\hfill \end{array}\right.$$

com $w_c$ constante e $w_0\left(t\right)$ uma função crescente com $w_0\left(t\right)\to w^{\infty }<\infty$ quando $t\to \infty$. Nesse problema a temperatura inicial $f\left(r\right)$ é a temperatura corporal que é constante e dada por $f_0$.

Uma característica interessante deste modelo é que ele é um sistema acoplado, isto é, as equações dependem uma da outra. Mas, pelo fato das constantes físicas do problema serem diferentes isso leva a um imbróglio na hora de encontrar as soluções. Uma estratégia comum utilizada para problemas como esse é desacoplar as equações utilizando a teoria de autovalores e autovetores na matriz que acopla o sistema, que neste caso é dada por

(6.2.3) $$\mathbf{P}=\left[\begin{array}{cc}-h_{tb}& h_{tb}\\ h_{tb}& -h_{tb}\end{array}\right].$$

É imediato perceber que essa matriz tem autovalores complexos, logo não é possível diagonalizá-la e realizar o procedimento de desacoplamento. Portanto, a estratégia de utilizar a transformada finita se faz bastanteapropriada nesta situação.

Nas seções a seguir iremos construir passo a passo a solução do problema definido pelas equações (6.2.1) sujeito as condições iniciais e de contorno (6.2.2). Começaremos encontrando os autovalores e autofunções do problema de Sturm--Liouville associado ao sistema de equações desejado. Em seguida, será definida a transformada finita para o problema (6.2.1)--(6.2.2). Finalmente, iremos aplicar a transformada finita definida nas equações do sistema para obter a solução do problema de valor inicial e de contorno dado.

6.2.1 Autovalores e autofunções

Considerando o problema com simetria radial, então o operador diferencial é dado por

(6.2.4) $${\Delta }_e=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}\right),$$

isto é, as funções $P\left(r\right),Q\left(r\right)$ e $W\left(r\right)$ são dadas por $P\left(r\right)=r^2$, $Q\left(r\right)=0$ e $W\left(r\right)=r^2$. Nesse problema o operador está associado a seguinte equação diferencial

(6.2.5) $${\left[r^2{R}^{\prime }\left(r\right)\right]}^{\prime }=-{\lambda }^2{r}^{2}R\left(r\right),$$

lembrando que neste caso o problema de Sturm--Liouville também é regular, pois $r\in (a,b)$ com $a>0$ Essa equação pode ser resolvida através da mudança de variáveis $R\left(r\right)=V\left(r\right)/r$, que a reduz a uma equação de segunda ordem linear e com coeficientes constantes cuja solução geral é dada na forma

(6.2.6) $$R\left(r\right)=c_1\frac{\text{sen}\left(\lambda r\right)}{r}+c_2\frac{cos\left(\lambda r\right)}{r},$$

observe que as soluções linearmente independentes são as funções de Bessel esféricas de ordem zero, $j_0\left(\lambda r\right)$ e $y_0\left(\lambda r\right)$, respectivamente ^[27]. Como esse problema tem condições de contorno baseadas no valor da temperatura nas extremidades, então iremos considerar que as condições de contorno do problema de Sturm--Liouville associado são homogêneas e dadas por $R\left(a\right)=R\left(b\right)=0.$ Nesses valores, temos o seguinte sistema

(6.2.7) $$\begin{array}{l}{c}_{1}sen\left(\lambda a\right)+c_{2}cos\left(\lambda a\right)=0\\ c_{1}sen\left(\lambda b\right)+c_{2}cos\left(\lambda b\right)=0.\end{array}$$

Além disso, considerando $cos\left(\lambda a\right)\ne 0$ então podemos escrever a constante $c_2$ em função da constante $c_1$ na forma

(6.2.8) $$c_2=-c_1\text{tg}\left(\lambda a\right).$$

Substituindo a equação (6.2.8) na equação (6.2.7), temos

(6.2.9) $$\frac{c_1}{cos\left(\lambda a\right)}\text{sen}\left(\lambda (b-a)\right)=0.$$

que nos dará a equação dos autovalores do problema, pois $c_1\ne 0$. Neste caso é fácil ver que os autovalores serão dados exatamente por

(6.2.10) $${\lambda }_m=\frac{m\pi }{b-a},$$

em que $m\in \mathbb{N}$. Consequentemente, temos que as autofunções associadas aos autovalores encontrados são dadas por

(6.2.11) $$R_m\left(r\right)=\frac{1}{r}\text{sen}\left({\lambda }_m\left(r-a\right)\right).$$

Diferentemente do exemplo na seção (6.1), perceba que neste exemplo o autovalor nulo não faz parte do conjunto de autovalores. Pois, o número ${\lambda }_0=0$ está associado a função $R_0\left(r\right)=0$.

6.2.2. Definindo a transformada

A transformada para esse problema será definida por

(6.2.12) $${\mathcal{H}}_m\left[g\right]={\int }_a^b{r}^{2}g\left(r\right)R_m\left(r\right)dr$$

em que a função $R_m\left(r\right)$ é dada pela equação (6.2.11).

De acordo com a equação (20), no operador diferencial ${\Delta }_e$ a transformada se comporta como

(6.2.13) $$\begin{array}{lll}{\mathcal{H}}_m\left[{\Delta }_{e}g\right]\hfill & =\hfill & {\lambda }_m^2{\mathcal{H}}_m\left[g\right]-{\left[r^{2}g\left(r\right)\frac{d{R}_m}{dr}\right]}_a^b\hfill \\ \hfill & =\hfill & -{\lambda }_m^2{\mathcal{H}}_m\left[g\right]-{\lambda }_m\left[{\left(-1\right)}^{m}bg\left(b\right)-ag\left(a\right)\right].\hfill \end{array}$$

Antes de aplicarmos essa transformada nas equações em (6.2.1) é interessante calcularmos o valor da transformada de algumas funções, em separado. Começamos com $g\left(r\right)=1$, temos que

(6.2.14) $$\begin{array}{lll}{\mathcal{H}}_m\left[1\right]\hfill & =\hfill & {\int }_a^b{r}^2\cdot 1\cdot \frac{sin\left({\lambda }_m\left(r-a\right)\right)}{r}dr\hfill \\ \hfill & =\hfill & {\int }_a^{b}rsin\left({\lambda }_m\left(r-a\right)\right)dr\hfill \\ \hfill & =\hfill & -\frac{1}{{\lambda }_m}\left[b{\left(-1\right)}^m-a\right].\hfill \end{array}$$

Além disso, para $g\left(r\right)=1/r^4$ temos

(6.2.15) $$\begin{array}{lll}\mathcal{H}\left[\frac{1}{r^4}\right]\hfill & =\hfill & {\int }_a^b{r}^2\cdot \frac{1}{r^4}\cdot \frac{sin\left({\lambda }_m\left(r-a\right)\right)}{r}dr\hfill \\ \hfill & =\hfill & {\int }_a^b\frac{sin\left({\lambda }_m\left(r-a\right)\right)}{r^3}dr\hfill \\ \hfill & =\hfill & \frac{1}{2}{\lambda }_m^{2}sen\left({\lambda }_{m}a\right)\left[Ci\left({\lambda }_{m}b\right)-Ci\left({\lambda }_{m}a\right)\right]\hfill \\ \hfill & -\hfill & \frac{1}{2}{\lambda }_m^{2}cos\left({\lambda }_{m}a\right)\left[Si\left({\lambda }_{m}b\right)-Si\left({\lambda }_{m}a\right)\right]\hfill \\ \hfill & -\hfill & \frac{1}{2}{\lambda }_m\left(\frac{{\left(-1\right)}^m}{b}-\frac{1}{a}\right),\hfill \end{array}$$

em que $\text{Si}\left(z\right)={\int }_0^z\frac{\text{sen}\left(x\right)}{x}dx$ é o Seno Integral e $\text{Ci}\left(z\right)=\gamma +ln\left(z\right)+{\int }_0^z\frac{cos\left(x\right)-1}{x}dx$ é o Cosseno Integral, com $\gamma =-{\int }_0^{\infty }{e}^{-x}ln\left(x\right)dx=0.57721\cdots$ a constante de Euler-Mascheroni ^[27].

Além disso, na solução final também será necessário o valor do produto interno que é dado por

(6.2.16) $$\begin{array}{lll}\langle R_m,R_m\rangle \hfill & =\hfill & {\int }_a^b{r}^2{\left[\frac{sen\left({\lambda }_m\left(r-a\right)\right)}{r^2}\right]}^{2}dr\hfill \\ \hfill & =\hfill & {\int }_a^{b}se{n}^2\left({\lambda }_m\left(r-a\right)\right)dr\hfill \\ \hfill & =\hfill & \frac{b-a}{2}.\hfill \end{array}$$

6.2.3. Resolvendo o sistema

Antes de aplicarmos a transformada definida pelas equações (6.2.12)--(6.2.13), primeiramente reescrevemos o sistema dado pela equação (6.2.1) na forma

(6.2.17) $$\begin{array}{l}\frac{\partial w_t}{\partial t}={\alpha }_t{\Delta }_e{w}_t-{\beta }_t\left(w_t-w_b\right)+n_t+\frac{{\theta }_t}{r^4}\\ \frac{\partial w_b}{\partial t}={\alpha }_b{\Delta }_e{w}_b-{\beta }_b\left(w_t-w_b\right)+n_b+\frac{{\theta }_b}{r^4},\end{array}$$

em que ${\alpha }_t=k_t/{\rho }_t{c}_t$,${\alpha }_b=k_b/{\rho }_b{c}_b$, ${\beta }_t=h_{tb}/{\rho }_t{c}_t$, ${\beta }_b=h_{tb}/{\rho }_b{c}_b$, ${\eta }_t=q_t/{\rho }_t{c}_t$, ${\eta }_b=q_b/{\rho }_b{c}_b$, ${\theta }_t={\psi }_t/{\rho }_t{c}_t$ e ${\theta }_b={\psi }_b/{\rho }_b{c}_b$. Assim, aplicando a transformada nas equações diferenciais em (6.2.17) temos que

(6.2.18) $$\begin{array}{lll}\frac{\partial {\tilde{w}}_t}{\partial t}\hfill & =\hfill & -{\alpha }_t{\lambda }_m^2{\tilde{w}}_t-{\beta }_t\left(\tilde{w}{}_t-{\tilde{w}}_b\right)-{\lambda }_{m}a{w}_0\left(t\right)\hfill \\ \hfill & +\hfill & {\left(-1\right)}^{m+1}b{w}_c+n_t{\mathcal{H}}_m\left[1\right]+{\theta }_t{\mathcal{H}}_m\left[\frac{1}{r^4}\right]\hfill \\ \frac{\partial {\tilde{w}}_b}{\partial t}\hfill & =\hfill & -{\alpha }_b{\lambda }_m^2{\tilde{w}}_b+{\beta }_b\left({\tilde{w}}_t-{\tilde{w}}_b\right)-{\lambda }_{m}a{w}_0\left(t\right)\hfill \\ \hfill & +\hfill & {\left(-1\right)}^{m-1}b{w}_c+\eta b{\mathcal{H}}_m\left[1\right]+{\theta }_b{\mathcal{H}}_m\left[\frac{1}{r^4}\right],\hfill \end{array}$$

para ${\tilde{w}}_t={\mathcal{H}}_m\left[w_t\right]$ e ${\tilde{w}}_b={\mathcal{H}}_m\left[w_b\right]$. Reescrevendo as equações em (6.2.18) em uma forma mais compacta temos

(6.2.19) $$\begin{array}{l}\frac{\partial {\tilde{w}}_t}{\partial t}+\left({\alpha }_t{\lambda }_m^2+{\beta }_t\right){\tilde{w}}_t-{\beta }_t{\tilde{w}}_b=a{\lambda }_m{w}_0\left(t\right)+g_t\left({\lambda }_m\right)\\ \frac{\partial {\tilde{w}}_b}{\partial t}+\left({\alpha }_b{\lambda }_m^2+{\beta }_b\right){\tilde{w}}_b-{\beta }_b{\tilde{w}}_t=a{\lambda }_m{w}_0\left(t\right)+g_b\left({\lambda }_m\right),\end{array}$$

em que $g_t\left({\lambda }_m\right)={\left(-1\right)}^{m+1}{\lambda }_{m}b{w}_c+\delta q_t{\mathcal{H}}_m\left[1\right]+\delta {\psi }_t{\mathcal{H}}_m\left[\frac{1}{r^4}\right]$ e $g_b\left({\lambda }_m\right)={\left(-1\right)}^{m+1}{\lambda }_{m}b{w}_c+\varepsilon q_b{\mathcal{H}}_m\left[1\right]+\varepsilon {\psi }_b{\mathcal{H}}_m\left[\frac{1}{r^4}\right]$ não dependem de $t$.

Perceba que as equações em (6.2.19) fornecem um sistema de equações diferenciais ordinárias acopladas na variável temporal. Assim, para resolver esse problema, podemos utilizar a transformada de Laplace que é definida na forma

(6.2.20) $$\mathcal{L}\left[g\right]={\int }_0^{\infty }g\left(t\right)e^{-st}ds.$$

É importante ressaltar que a utilização da transformada de Laplace facilita a obtenção da solução final, porém não é a única alternativa para a solução do problema transiente. É possível realizar um desacoplamento do sistema de equações ordinárias dado pelas equações em (6.2.19) utilizando autovalores e autovetores para obter a mesma solução ao final.

Logo, aplicando a transformada de Laplace em cada uma das equações em (6.2.19) se obtém

(6.2.21) $$\begin{array}{ccc}& & \left(s+{\alpha }_t{\lambda }_m^2+{\beta }_t\right){\widehat{w}}_t-{\beta }_t{\widehat{w}}_b=a{\lambda }_m{\widehat{w}}_0\left(s\right)\hfill \\ & & +\frac{1}{s}{g}_t\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right)\hfill \\ & & \left(s+{\alpha }_b{\lambda }_m^2+{\beta }_b\right){\widehat{w}}_b-{\beta }_b{\widehat{w}}_t=a{\lambda }_m{\widehat{w}}_0\left(s\right)\hfill \\ & & +\frac{1}{s}{g}_b\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right),\hfill \end{array}$$

em que ${\widehat{w}}_t=\mathcal{L}\left[{\tilde{w}}_t\right]$, ${\widehat{w}}_b=\mathcal{L}\left[{\tilde{w}}_b\right]$ e $\tilde{f}\left({\lambda }_m\right)={\mathcal{H}}_m\left[f\right]$.

Para resolver o sistema dado pelas equações em (6.2.21), vamos escrevê-lo em notação matricial na forma

(6.2.22) $$\mathbf{M}\mathbf{v}=\mathbf{z},$$

em que a matriz dos coeficientes é dada por

(6.2.23) $$\mathbf{M}=\left[\begin{array}{cc}s+{\alpha }_t{\lambda }_m^2+{\beta }_t& -{\beta }_t\\ -{\beta }_b& s+{\alpha }_b{\lambda }_m^2+{\beta }_b\end{array}\right],$$

e os vetores dados por $\mathbf{v}={\left[{\widehat{w}}_t,{\widehat{w}}_b\right]}^T$ e

(6.2.24) $$\mathbf{z}=\left[\begin{array}{c}a{\lambda }_m{\widehat{w}}_0\left(s\right)+\frac{1}{s}{g}_t\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right)\\ a{\lambda }_m{\widehat{w}}_0\left(s\right)+\frac{1}{s}{g}_b\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right)\end{array}\right].$$

A matriz não-singular $\mathbf{M}$ tem sua inversa dada por

(6.2.25) $${\mathbf{M}}^{-1=}\frac{1}{det\mathbf{M}}\left[\begin{array}{cc}s+{\alpha }_b{\lambda }_m^2+{\beta }_b& {\beta }_t\\ {\beta }_b& s+{\alpha }_t{\lambda }_m^2+{\beta }_t\end{array}\right]$$

e seu determinante, det M, pode ser escrito como

(6.2.26) $$\begin{array}{ccc}\hfill det\mathbf{M}& =& \left[s+({\lambda }_m^2{\alpha }_t+{\beta }_t)\right]\left[s+({\lambda }_m^2{\alpha }_b+{\beta }_b)\right]-{\beta }_b{\beta }_t\hfill \\ & =& s^2+s\left({\lambda }_m^2{\alpha }_b+{\beta }_b+{\lambda }_m^2{\alpha }_t+{\beta }_t\right)\hfill \\ & +& ({\lambda }_m^2{\alpha }_t+{\beta }_t)({\lambda }_m^2{\alpha }_b+{\beta }_b)-{\beta }_b{\beta }_t\hfill \\ & =& {\left[s+\frac{1}{2}({\lambda }_m^2{\alpha }_b+{\beta }_b+{\lambda }_m^2{\alpha }_t+{\beta }_t)\right]}^2\hfill \\ & -& \frac{1}{4}\left\{\left[({\lambda }_m^2{\alpha }_b+{\beta }_b)-({\lambda }_m^2{\alpha }_t+{\beta }_t)\right]{}^2+4{\beta }_b{\beta }_t\right\}\hfill \\ & =& {(s+{\psi }_m)}^2-{\phi }_m^2,\hfill \end{array}$$

em que ${\psi }_m=\frac{1}{2}({\psi }_{m,b}+{\psi }_{m,t})$, ${\phi }_m^2=\frac{1}{4}\left[\left({\psi }_{m,b}-{\psi }_{m,t}\right){}^2+4{\beta }_b{\beta }_t\right]$, ${\psi }_{m,b}={\lambda }_m^2{\alpha }_b+{\beta }_b$ e ${\psi }_{m,t}={\lambda }_m^2{\alpha }_t+{\beta }_t$. Então, de forma compacta, nós temos que a solução do sistema linear é dado por

(6.2.27) $$\begin{array}{l}{\stackrel{\circ }{w}}_t=\frac{1}{det\text{M}}\left[\frac{A_m}{s}+B_{m}s+C_m+D_m\left(s\right)s+E_m\left(s\right)\right]\\ {\stackrel{\circ }{w}}_b=\frac{1}{det\text{M}}\left[\frac{F_m}{s}+G_{m}s+H_m+I_m\left(s\right)s+J_m\left(s\right)\right],\end{array}$$

em que os coeficientes são

(6.2.28) $$\begin{array}{ccc}\hfill A_m& =& {\psi }_b{g}_t\left({\lambda }_m\right)+{\beta }_t{g}_b\left({\lambda }_m\right)\hfill \end{array}$$

(6.2.29) $$\begin{array}{ccc}\hfill B_m& =& \tilde{f}\left({\lambda }_m\right)\hfill \end{array}$$

(6.2.30) $$\begin{array}{ccc}\hfill C_m& =& g_t\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right)\left({\psi }_b+{\beta }_t\right)\hfill \end{array}$$

(6.2.31) $$\begin{array}{ccc}\hfill D_m\left(s\right)& =& a{\lambda }_m{\widehat{w}}_0\left(s\right)\hfill \end{array}$$

(6.2.32) $$\begin{array}{ccc}\hfill E_m\left(s\right)& =& a{\lambda }_m{\widehat{w}}_0\left(s\right)\left({\psi }_b+{\beta }_t\right)\hfill \end{array}$$

(6.2.33) $$\begin{array}{ccc}\hfill F_m& =& {\psi }_t{g}_b\left({\lambda }_m\right)+{\beta }_b{g}_t\left({\lambda }_m\right)\hfill \end{array}$$

(6.2.34) $$\begin{array}{ccc}\hfill G_m& =& \tilde{f}\left({\lambda }_m\right)\hfill \end{array}$$

(6.2.35) $$\begin{array}{ccc}\hfill H_m& =& g_b\left({\lambda }_m\right)+\tilde{f}\left({\lambda }_m\right)\left({\psi }_t+{\beta }_b\right)\hfill \end{array}$$

(6.2.36) $$\begin{array}{ccc}\hfill I_m\left(s\right)& =& a{\lambda }_m{\widehat{w}}_0\left(s\right)\hfill \end{array}$$

(6.2.37) $$\begin{array}{ccc}\hfill J_m\left(s\right)& =& a{\lambda }_m{\widehat{w}}_0\left(s\right)\left({\psi }_t+{\beta }_b\right).\hfill \end{array}$$

Considerando as seguintes transformadas inversas de Laplace

(6.2.38) $$\begin{array}{c}\hfill {\displaystyle e^{-{\psi }_{m}t}cosh\left({\phi }_{m}t\right)}={\mathcal{L}}^{-1}\left[\frac{s+{\psi }_m}{{(s+{\psi }_m)}^2-{\phi }_m^2}\right]\\ \hfill {\displaystyle e^{-{\psi }_{m}t}sinh\left({\phi }_{m}t\right)}={\mathcal{L}}^{-1}\left[\frac{{\phi }_m}{{(s+{\psi }_m)}^2-{\phi }_m^2}\right]\end{array}$$

temos então as seguintes soluções do sistema de equações diferenciais (6.2.1) associado as condições (6.2.2) que convergem no domínio $(r,t)\in (a,b)\times [0,\infty )$ e formalmente são dadas por

(6.2.39) $$\begin{array}{ccc}\hfill w_t(r,t)& =& \sum _{m=1}^{\infty }\left\{\frac{A_m}{{\psi }_m^2-{\phi }_m^2}+\left(B_m-\frac{A_m}{{\psi }_m^2-{\phi }_m^2}\right)e^{-{\psi }_{m}t}cosh\left({\phi }_{m}t\right)\right.\hfill \\ & +& \left.\frac{1}{{\phi }_m}\left(C_m-{\psi }_m{B}_m-\frac{{\psi }_m{A}_m}{{\psi }_m^2-{\phi }_m^2}\right)e^{-{\psi }_{m}t}sinh\left({\phi }_{m}t\right)\right\}{N}_m\left(r\right)\hfill \\ & +& \sum _{m=1}^{\infty }\left\{{\int }_0^t{D}_m(t-\tau )e^{-{\psi }_m\tau }cosh\left({\phi }_m\tau \right)d\tau \right.\hfill \\ & +& \left.{\int }_0^t\frac{1}{{\phi }_m}\left[E_m(t-\tau )-{\psi }_m{D}_m(t-\tau )\right]e^{-{\psi }_m\tau }sinh\left({\phi }_m\tau \right)d\tau \right\}{N}_m\left(r\right)\hfill \end{array}$$

(6.2.40) $$\begin{array}{ccc}\hfill w_b(r,t)& =& \sum _{m=1}^{\infty }\left\{\frac{F_m}{{\psi }_m^2-{\phi }_m^2}+\left(G_m-\frac{F_m}{{\psi }_m^2-{\phi }_m^2}\right)e^{-{\psi }_{m}t}cosh\left({\phi }_{m}t\right)\right.\hfill \\ & +& \left.\frac{1}{{\phi }_m}\left(H_m-{\psi }_m{G}_m-\frac{{\psi }_m{F}_m}{{\psi }_m^2-{\phi }_m^2}\right)e^{-{\psi }_{m}t}sinh\left({\phi }_{m}t\right)\right\}{N}_m\left(r\right)\hfill \\ & +& \sum _{m=1}^{\infty }\left\{{\int }_0^t{I}_m(t-\tau )e^{-{\psi }_m\tau }cosh\left({\phi }_m\tau \right)d\tau \right.\hfill \\ & +& \left.{\int }_0^t\frac{1}{{\phi }_m}\left[J_m(t-\tau )-{\psi }_m{I}_m(t-\tau )\right]e^{-{\psi }_m\tau }sinh\left({\phi }_m\tau \right)d\tau \right\}{N}_m\left(r\right),\hfill \end{array}$$

em que as autofunções normalizadas são dadas por

(6.2.41) $$N_m\left(r\right)=\frac{2}{b-a}\text{sen}\left(\frac{m\pi }{b-a}\left(r-a\right)\right)$$

e os coeficientes pelas equações (6.2.28)--(6.2.37).

É interessante observar nessa solução que ${\psi }_m>{\phi }_m$, caso contrário as soluções $w_t,w_b\to \infty$ quando $t\to \infty$. Além disso, é possível resolver esse problema de uma forma distinta decompondo a temperatura de cada fase como a soma de uma solução permanente e uma transiente. Essa segunda forma de resolver o problema é pertinente, pois as funções dadas pelas equações (6.2.39) e (6.2.40) satisfazem $w_t(a,t)=w_b(a,t)=0$ e $w_t(b,t)=w_b(b,t)=0$, isto é, não convergem para as condições de contorno nos extremos do intervalo $[a,b]$. Enquanto com a decomposição sugerida, entre solução permanente e solução transiente, a solução final da equação não sofrerá o problema de convergência nos extremos e convergirá em todo domínio $(r,t)\in [a,b]\times [0,\infty )$.

7. Considerações finais e perspectivas

Solucionar problemas relacionados as equações diferenciais é um trabalho que faz parte do cotidiano de físicos, matemáticos e engenheiros. Dessa forma, o acesso a diferentes técnicas de solução é extremamente importante para o trabalho desses profissionais. Nesse sentido, foi apresentada neste trabalho a técnica de solução de problemas de valores inicial e de contorno por meio da construção de transformações integrais finitas baseadas nas soluções de problemas de Sturm--Liouville. Este tipo de método é interessante, pois a solução é obtida de uma forma construtiva. Ele proporciona uma visão distinta sobre a expansão por autofunções levando a uma forma mais intuitiva e operacional de se resolver os problemas, semelhante as regras operacionais das Transformadas de Laplace e Fourier. Além de apresentar uma metodologia diferente da usual função de Green.

Assim, no primeiro momento, foi apresentada uma breve revisão acerca da equação do calor, que foi a base para a descrição do método, e dos problemas de Sturm--Liouville regulares. Na sequência, foi apresentado como estabelecer uma transformação integral apropriada ao problema que se deseja solucionar. Por fim, foram discutidas as soluções de dois problemas distintos relacionados a pesquisas que estão sendo desenvolvidas na Universidade de Brasília. As soluções foram construídas de forma didática, primeiramente para mostrar a importância desta técnica, e também com a perspectiva de que o leitor seja capaz de reproduzir e procurar soluções de outros problemas de valor inicial e de contorno utilizando a transformada finita.

Concluímos dizendo que a técnica aqui mostrada e exemplificada não tem como objetivo se apresentar como uma técnica definitiva para solução de problemas de valor inicial e de contorno para equações diferenciais parciais lineares. Mas, sim, como mais uma ferramenta alternativa que pode ser adicionada a coleção de técnicas que o leitor possui para a solução de problemas nessa forma e também como uma outra visão, pouco explorada nos livros texto, sobre a já conhecida técnica da expansão por autofunções.

Referências

• ^[1]
  F. Cajori. Uma História da Matemática (Ciências Moderna, Rio de Janeiro, 2007).
• ^[2]
  J. Stillwell, Mathematics and Its History (Springer, New York, 2010).
• ^[3]
  W. E. Boyce e R.C. DiPrima, Equações Diferenciais Elementares e Problemas de Valores de Contorno (LTC, Rio de Janeiro, 1979).
• ^[4]
  R. I. Junior, V. M. Iório, Equações Diferenciais Parciais: uma Introdução (IMPA, Rio de Janeiro, 2013).
• ^[5]
  J. B. Filho, E. C. Oliveira and H. G. Pavão. Revista Brasileira de Ensino de Física, 6, 2 (1984).
• ^[6]
  E. C. Oliveira and M. Tygel. Métodos Matemáticos para Engenharia (SBM, Rio de Janeiro, 2010)
• ^[7]
  H. K. Brown. Resolution of boundary value problems by means of the finite Fourier transformation; general vibration of a string. J. Appl. Phys. 14, 609 (1944).
• ^[8]
  I. Roettinger. A generalization of the finite Fourier transformation and applications. Quart. Appl. Math. 5, 298 (1947).
• ^[9]
  I. N. Sneddon. III. Finite Hankel Transforms. Lond. Edinb. Dubl. Phil. Mag. 37, 17 (1946).
• ^[10]
  G. Cinelli. An extension of the finite Hankel transform and applications. Int. J. Eng. Sci. 3, 539 (1965).
• ^[11]
  H. S. Dunn. A generalization of the Laplace transform. Math. Proc. Camb. Philos. Soc. 63, 155 (1967).
• ^[12]
  R. Datko. Applications of the Finite Laplace Transform to Linear Control Problems. SIAM J. Control Optim. 18, 386 (1980).
• ^[13]
  G. B. Arfken, H. J. Weber. Mathematical Methods for Physicists (Elsevier, Oxford, 2005).
• ^[14]
  L. Debnath and D. Bhatta. Integral Transforms and Their Applications (Chapman & Hall, Boca Raton, 2007).
• ^[15]
  J. Irving, Mullineux, N. Mathematics in Physics and Engineering (Academic Press, New York, 1959).
• ^[16]
  M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, New York, 1972).
• ^[17]
  M. A. Al-Gwaiz. Sturm-Liouville Theory and Its Applications (Springer, London, 2008).
• ^[18]
  J. Muscat. Functional Analisys (Springer, Switzerland, 2014).
• ^[19]
  A. P. C. P. Nunes, G. P. Milhomem, V. C. Rispoli and A. Andrianov, in Proceedings of the XXXVIII Ibero-Latin American Congress on Computational Methods in Engineering, Florianópolis, 2017, editado por P.O. Faria, R.H. Lopez, L.F.F. Miguel, W.J.S. Gomes, M. Noronha (ABMEC, Florianópolis, 2017), p. 1.
• ^[20]
  M. Monteiro, S. F. R. Rosa, B. Motta, G. Anjos, M. P. Marques. The Use of Radiofrequency for Hepatocellular Carcinoma Ablation: an Update Review and Perspectives. Int. J. Biosen. Bioelectron. 1, 55 (2017).
• ^[21]
  T. Peng, D. P. O'Neill, S. J. Payne. A two-equation coupled system for determination of liver tissue temperature during thermal ablation. Int. J. Heat Mass Transf. 54, 2100 (2010).
• ^[22]
  Y. Pinchover and J. Rubinstein. An Introduction to Partial Differential Equations. (Cambridge University Press, Cambridge, 2005).
• ^[23]
  T. J. Ypma. Historical Development of the Newton--Raphson Method. SIAM Review 37, 531 (1995).
• ^[24]
  F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions (Cambridge University Press, Cambridge, 2010).
• ^[25]
  B. Zhang, M. A. J. Moser, E. M. Zhang, Y. Luo and W. Zhang, in Proceedings of the 51st IFMBE World Congress on Medical Physics and Biomedical Engineering, (Toronto, Canadá, 2015).
• ^[26]
  K. Wang, F. Tavakkoli, S. Wang and K. Vafai. Analysis and analytical characterization of bioheat transfer during radiofrequency ablation. J. Biomech. 48, 930 (2015).
• ^[27]
  I. N. Bronshtein, K. A. Semendyayev, G. Musiol and H. Mühlig. Handbook of Mathematics (Springer-Verlag, Berlin, 2015).

Datas de Publicação

• Publicação nesta coleção
  2018

Histórico

• Recebido
  14 Dez 2017
• Revisado
  19 Mar 2018
• Aceito
  03 Abr 2018