sujeita a uma condição inicial c(x, 0)=1 e uma reação de depleção com cinética de primeira ordem em cada contorno. Ou seja,
\boldsymbol{n} \cdot \mathrm{J}=-k c
Onde \boldsymbol{n} é um vetor normal que aponta para fora do domínio, enquanto:
\mathbf{J}=-D \boldsymbol{\nabla} c
é o vetor fluxo difusivo. Determine o sistema de equações diferenciais ordinárias que você precisaria resolver para usar o Método das Linhas e Diferenças Finitas Centradas.
Passo 1
O que o enunciado dá e o que pede, em linguagem clara e lista de dados úteis:
Condição inicial: c(x,0)=1 \quad\text{para todo }x\text{ no domínio.}
Condição de contorno nas duas extremidades (condição de fluxo de superfície do tipo Robin): \boldsymbol{n}\cdot\mathbf{J}=-k\,c,
com \mathbf{J}=-D\boldsymbol{\nabla}c. Em 1D isso será convertido em uma relação entre \partial c/\partial x e c nas extremidades.
Não foi dado explicitamente o comprimento do domínio; para ser concreto vou considerar o domínio unidimensional [0,L]. Isso é padrão em problemas assim; o procedimento não muda muito se L for dado.
Objetivo pedido: discretizar espacialmente por diferenças finitas centradas (Método das Linhas) e escrever o sistema de equações ordinárias em tempo que deve ser resolvido (i.e., obter \dot{c}_i(t) para cada nó espacial).
Passo 2
Escolha de malha espacial e notação (para depois escrever as fórmulas de diferenças):
Dividimos o intervalo [0,L] em N subintervalos de tamanho \Delta x=L/N.
Definimos os nós \,x_i = i\,\Delta x\,,\; i=0,1,\dots,N.
Denotamos a aproximação da concentração no nó i por \,c_i(t)\approx c(x_i,t).
O Método das Linhas: discretizamos as derivadas espaciais (aqui com diferenças finitas centradas) mantendo a derivada temporal contínua. Assim obtemos um sistema de ODEs para as incógnitas \{c_i(t)\}_{i=0}^{N}.
Fórmula central padrão para segunda derivada (em um nó interior): \frac{\partial^2 c}{\partial x^2}\Big|_{x_i}\approx \frac{c_{i-1}-2c_i+c_{i+1}}{\Delta x^{2}}.
Vou explicar e usar essa fórmula para os nós interiores; nas fronteiras precisamos usar a condição de contorno para eliminar pontos “fantasmas” se usarmos formula central também na fronteira.
Passo 3
Equações para os nós interiores (não nas extremidades):
Para nós interiores i=1,2,\dots,N-1 aplicamos diretamente a aproximação central para \partial^2 c/\partial x^2. Substituindo na PDE obtemos as ODEs espaciais: \frac{d c_i}{d t}=D\frac{c_{i-1}-2c_i+c_{i+1}}{\Delta x^{2}},\qquad i=1,\dots,N-1.
Aqui cada variável:
\,c_i(t) é a concentração aproximada no nó i;
D é o coeficiente de difusão (constante);
\Delta x é o espaçamento espacial.
Interpretação prática: cada \dot c_i é uma combinação linear das concentrações vizinhas. Isso é exatamente o passo do Método das Linhas: transformamos a parte espacial numa operação matricial sobre o vetor de c's.
Passo 4
Tradução da condição de contorno para 1D e expressão para as derivadas em \,x=0 e \,x=L:
Em 1D o fluxo é J_x=-D\,\partial c/\partial x.
No nó esquerdo x=0 o vetor normal sai do domínio para a esquerda, então n=-1. A condição \boldsymbol{n}\cdot\mathbf{J}=-k c dá -J_x|_{0}=-k c(0). Como J_x=-D\partial c/\partial x, isso resulta em:\frac{\partial c}{\partial x}\Big|_{x=0}=-\frac{k}{D}\,c(0).
No nó direito x=L o vetor normal aponta para a direita, n=+1. A condição dá J_x|_{L}=-k c(L). Novamente usando J_x=-D\partial c/\partial x:\frac{\partial c}{\partial x}\Big|_{x=L}=+\frac{k}{D}\,c(L).
Essas duas fórmulas são condições de contorno do tipo Robin (mistura de derivada e valor).
Passo 5
Como aplicar diferenças centradas nas fronteiras usando pontos fantasmas e eliminar os fantasmas com as condições de contorno:
Para manter a diferença central consistente também na segunda derivada nos nós i=0 e i=N, introduzimos nós fantasmas c_{-1} (à esquerda de x_0) e c_{N+1} (à direita de x_N). A fórmula central para a derivada espacial em fronteira é, por exemplo, para x=0:\frac{\partial c}{\partial x}\Big|_{x=0}\approx\frac{c_{1}-c_{-1}}{2\Delta x}.Igualando à condição de contorno \partial c/\partial x|_{0}=-\dfrac{k}{D}c_0 obtemos uma expressão para o ponto fantasma:\frac{c_{1}-c_{-1}}{2\Delta x}=-\frac{k}{D}c_0 \quad\Longrightarrow\quad c_{-1}=c_{1}+2\Delta x\frac{k}{D}c_0.
Analogamente, na fronteira direita usamos\frac{\partial c}{\partial x}\Big|_{x=L}\approx\frac{c_{N+1}-c_{N-1}}{2\Delta x}=\frac{k}{D}c_Nque dác_{N+1}=c_{N-1}+2\Delta x\frac{k}{D}c_N.
Agora podemos usar essas expressões para eliminar c_{-1} e c_{N+1} nas fórmulas de segunda derivada central nos nós de contorno.
Passo 6
Derivação explícita das ODEs para os nós de contorno substituindo os pontos fantasmas:
Para i=0 (nó esquerdo) a segunda derivada central é \dfrac{c_{-1}-2c_0+c_{1}}{\Delta x^{2}}. Substituindo c_{-1}=c_{1}+2\Delta x\frac{k}{D}c_0 obtemos:\frac{c_{-1}-2c_0+c_1}{\Delta x^2} =\frac{(c_1+2\Delta x\frac{k}{D}c_0)-2c_0+c_1}{\Delta x^2} =\frac{2(c_1-c_0)}{\Delta x^2}+\frac{2k}{D\Delta x}c_0.Multiplicando por D (pois \partial c/\partial t = D \partial^2 c/\partial x^2) a ODE no nó 0 fica:\frac{d c_0}{d t} =\frac{2D}{\Delta x^2}\,(c_1-c_0)+\frac{2k}{\Delta x}\,c_0.
Para i=N (nó direito) a segunda derivada central é \dfrac{c_{N-1}-2c_N+c_{N+1}}{\Delta x^{2}}. Substituindo c_{N+1}=c_{N-1}+2\Delta x\frac{k}{D}c_N obtemos:\frac{c_{N-1}-2c_N+c_{N+1}}{\Delta x^2} =\frac{2(c_{N-1}-c_N)}{\Delta x^2}+\frac{2k}{D\Delta x}c_N,e a ODE no nó N é:\frac{d c_N}{d t} =\frac{2D}{\Delta x^2}\,(c_{N-1}-c_N)+\frac{2k}{\Delta x}\,c_N.
Observação rápida de interpretação: os termos com \dfrac{2D}{\Delta x^2} representam acoplamento com o vizinho interior; os termos com \dfrac{2k}{\Delta x} vêm da condição de fluxo proporcional a c e aparecem como contribuição local nas equações de contorno.
Passo 7
Resumo organizando todo o sistema de ODEs (pronto para ser integrado no tempo) e condição inicial:
Nós interiores i=1,\dots,N-1: \frac{d c_i}{d t}=D\frac{c_{i-1}-2c_i+c_{i+1}}{\Delta x^{2}},\qquad i=1,\dots,N-1.
Nó esquerdo i=0: \frac{d c_0}{d t}=\frac{2D}{\Delta x^{2}}\,(c_{1}-c_0)+\frac{2k}{\Delta x}\,c_0.
Nó direito i=N: \frac{d c_N}{d t}=\frac{2D}{\Delta x^{2}}\,(c_{N-1}-c_N)+\frac{2k}{\Delta x}\,c_N.
Em forma matricial compacta (vetor \mathbf{c}=[c_0,\dots,c_N]^\top) podemos escrever \dfrac{d\mathbf{c}}{dt}=A\,\mathbf{c} onde A é uma matriz tridiagonal modificada; os elementos são dados pelas fórmulas acima (termos diagonais e sub/sobrediagonais explícitos para os nós de contorno e valores habituais \dfrac{D}{\Delta x^2}\{-2,1,1\} para nós interiores).