Problema de transporte

Use o botão abaixo para reportar erros ou dar sugestões.

Computação Científica - Um Livro Colaborativo

13.2 Problema de transporte

Exemplo 13.2.1. Resolva o problema dado pela equação estacionária do transporte unidimensional com espalhamento isotrópico: $\begin{array}{l} μ \frac{\partial I}{\partial y} + λ I & = \frac{σ}{2} \int_{- 1}^{1} I (y, μ) d μ + Q (y, μ), & y \in (0, L), μ \in [- 1, 1], \\ I (0, μ) & = B_{0}, & μ > 0, \\ I (L, μ) & = B_{L}, & μ < 0 . \end{array}$

Nesse modelo, $I (y, μ)$ é uma intensidade radiativa, $I (y) = \frac{1}{2} \int_{- 1}^{1} I (y, μ) d μ$

é o fluxo escalar, $λ$ é o coeficiente de absorção, $σ$ é o coeficiente de espalhamento, $Q (y, μ)$ é uma fonte, $B_{0}$ e $B_{L}$ são os fluxos de entrada na fronteira, $y \in [0, L]$ é a variável espacial e $μ \in [- 1, 1]$ é o cosseno do ângulo entre o raio e o eixo $y$ .

Começamos discretizando à variável angular $μ$ com um método chamado de Ordenadas Discretas. Escolhemos a quadratura numérica de Gauss-Legendre com um número par de abscissas $μ_{i}$ e pesos $ω_{i}$ e discretizamos a equação de forma a transformar a integral em um somatório. Este procedimento transforma o problema em um sistema de EDO’s na variável $y$ : $\begin{array}{l} μ_{i} \frac{d I_{i}}{d y} + λ I_{i} & = \frac{σ}{2} \sum_{j = - M}^{M} ω_{j} I_{j} (y) + Q_{i}, & y \in (0, L), \\ I_{i} (0) & = B_{0}, & μ_{i} > 0, \\ I_{i} (L) & = B_{L}, & μ_{i} < 0, \end{array}$

onde $i = - M, . . ., M$ . Agora, separaramos o sistema de EDO’s em dois: um para $μ_{i} < 0$ e outro para $μ_{i} > 0$ . Portanto, $\begin{array}{l} - μ_{i} \frac{d I_{- i}}{d y} + λ I_{- i} & = \frac{σ}{2} \sum_{j = 1}^{M} ω_{j} (I_{- j} + I_{j}) + Q_{- i}, & i = 1, 2, . . ., M, \\ I_{- i} (L) & = B_{L}, & i = 1, 2, . . ., M, \end{array}$

e $\begin{array}{l} μ_{i} \frac{d I_{i}}{d y} + λ I_{i} & = \frac{σ}{2} \sum_{j = 1}^{M} ω_{j} (I_{- j} + I_{j}) + Q_{i}, & i = 1, 2, . . ., M, \\ I_{i} (0) & = B_{0}, & i = 1, 2, . . ., M . \end{array}$

Agora, vamos discretizar à variável espacial $y$ . Dividimos o intervalo $[0, L]$ em $N$ subintervalos de tamanho $h = \frac{L}{N}$ , e construímos a malha $y_{k} = (k - 1) h$ , $k = 1, 2, . . ., N + 1$ . Suponhamos que em um intervalo da malha $[y_{k}, y_{k + 1}]$ o lado direito das equações seja constante e conhecido. Vamos denotar a aproximação para o lado direito da equação de $S_{i}^{k + 1 ∕ 2}$ . Logo, neste intervalo temos: $\begin{array}{l} \frac{d I_{- i}}{d y} - \frac{λ}{μ_{i}} I_{- i} & = - \frac{1}{μ_{i}} S_{- i}^{k + 1 ∕ 2}, & i = 1, 2, . . ., M, & y_{k} \leq y \leq y_{k + 1}, \\ \frac{d I_{i}}{d y} + \frac{λ}{μ_{i}} I_{i} & = \frac{1}{μ_{i}} S_{i}^{k + 1 ∕ 2}, & i = 1, 2, . . ., M, & y_{k} \leq y \leq y_{k + 1} . \end{array}$

As duas equações podem ser resolvidas pelo método do fator integrante: $\begin{array}{l} e^{- \frac{λ}{μ_{i}} y_{k + 1}} I_{- i}^{k + 1} - e^{- \frac{λ}{μ_{i}} y_{k}} I_{- i}^{k} & = \frac{S_{- i}^{k + 1 ∕ 2}}{λ} (e^{- \frac{λ}{μ_{i}} y_{k + 1}} - e^{- \frac{λ}{μ_{i}} y_{k}}), & i = 1, 2, . . ., M, & k = 1, 2, . . ., N, \\ e^{\frac{λ}{μ_{i}} y_{k + 1}} I_{i}^{k + 1} - e^{\frac{λ}{μ_{i}} y_{k}} I_{i}^{k} & = \frac{S_{i}^{k + 1 ∕ 2}}{λ} (e^{\frac{λ}{μ_{i}} y_{k + 1}} - e^{\frac{λ}{μ_{i}} y_{k}}), & i = 1, 2, . . ., M, & k = 1, 2, . . ., N . \end{array}$

As condições de contorno nos dão: $\begin{array}{l} I_{- i}^{N + 1} & = B_{L}, & i = 1, 2, . . ., M, \\ I_{i}^{1} & = B_{0}, & i = 1, 2, . . ., M . \end{array}$

Logo, podemos obter $I_{- i}^{k}$ e $I_{i}^{k + 1}$ de maneira recursiva, ou seja, $\begin{array}{l} I_{- i}^{k} & = e^{\frac{λ}{μ_{i}} (y_{k} - y_{k + 1})} I_{- i}^{k + 1} + \frac{S_{- i}^{k + 1 ∕ 2}}{λ} (1 - e^{\frac{λ}{μ_{i}} (y_{k} - y_{k + 1})}), & i = 1, 2, . . ., M, & k = N, N - 1, . . ., 1, \\ I_{i}^{k + 1} & = e^{\frac{λ}{μ_{i}} (y_{k} - y_{k + 1})} I_{i}^{k} + \frac{S_{i}^{k + 1 ∕ 2}}{λ} (1 - e^{\frac{λ}{μ_{i}} (y_{k} - y_{k + 1})}), & i = 1, 2, . . ., M, & k = 1, 2, . . ., N . \end{array}$

Damos um chute inicial para a solução e calculamos $S_{i}^{k + 1 ∕ 2}$ com as seguintes médias: $\begin{array}{l} S_{\pm i}^{k + 1 ∕ 2} & = \frac{1}{2} [\frac{σ}{2} \sum_{j = 1}^{M} ω_{j} (I_{- j}^{k} + I_{j}^{k}) + Q_{\pm i}^{k} + \frac{σ}{2} \sum_{j = 1}^{M} ω_{j} (I_{- j}^{k + 1} + I_{j}^{k + 1}) + Q_{\pm i}^{k + 1}] . \end{array}$

A solução aproximada é obtida quando houver convergência do processo iterativo: Calculamos $S_{\pm i}^{k + 1 ∕ 2} = 0$ , depois $I_{- i}^{N + 1}$ , $I_{- i}^{N}$ , $I_{- i}^{N - 1}$ , $\dots$ , $I_{- i}^{1}$ e $I_{i}^{1}$ , $I_{i}^{2}$ , $\dots$ , $I_{i}^{N + 1}$ . Depois recalculamos $S_{\pm i}^{k + 1 ∕ 2}$ e $I_{\pm i}^{k}$ .

Para implementar esse algoritmo, fazemos algumas considerações:

É necessário entender o algoritmo antes de começar a programar. Embora seja uma questão óbvia, a ansiedade pode nos trair.
É importante pensar na estrutura do código antes de programar: tipo das principais variáveis, formas de armazenamento (ponteiros, matrizes), variáveis globais, estruturas, funções, pacotes, etc.
Não é indicado programar tudo de uma só vez. Faça por pedaços e valide cada função independentemente.
Compile o código seguidamente para não acumular erros de compilação.

Começamos implemetando uma função que calcula pesos e abscissas para a quadratura de Gauss-Legendre:

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#define M 10

double *x,*w;

double GL()
{
  size_t i,flag;
  gsl_integration_glfixed_table * table;
  table = gsl_integration_glfixed_table_alloc (M);
  for (i=0; i<M;i++)
  {
    flag = gsl_integration_glfixed_point (-1, 1, i, &x[i], &w[i], table);
    //printf("flag=%zu, %zu %f %f\n",flag, i, x[i], w[i]);
  }
  gsl_integration_glfixed_table_free(table);
}
double teste(double x)
{
  return x*x;
}
void main()
{
  int i;
  if (((x=malloc(sizeof(double)*M))==NULL)||((w=malloc(sizeof(double)*M))==NULL))
  {
  printf("Erro ao alocar memória");
  return;
  }
  GL(x,w);
  double result=0;
  for (i=0; i<M;i++) result+=w[i]*teste(x[i]);
  printf("O valor da integral é %f\n",result);
}

Depois de testar essa rotina integrando várias funções, vamos um passo adiante. Definimos duas variáveis globais para salvar pesos e abscissas e alocamos espaço para salvar a solução numa variável tipo ponteiro para ponteiro.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#define M 10
#define N 10

double *x,*w;

double GL()
{
  size_t i,flag;
  gsl_integration_glfixed_table * table;
  table = gsl_integration_glfixed_table_alloc (2*M);
  for (i=0; i<2*M;i++)
  {
    flag = gsl_integration_glfixed_point (-1, 1, i, &x[i], &w[i], table);
    //printf("flag=%zu, %zu %f %f\n",flag, i, x[i], w[i]);
  }
  gsl_integration_glfixed_table_free(table);
}

void main()
{
  int i;
  double **I;
  if (((x=malloc(sizeof(double)*2*M))==NULL)||((w=malloc(sizeof(double)*2*M))==NULL))
  {
  printf("Erro ao alocar memória");
  return;
  }
  if ((I = calloc(N+1,sizeof(double)))==NULL)
  {
  printf("Erro ao alocar memória");
  return;
  }
  for (i=0;i<N+1;i++) if ((I[i] = calloc(2*M,sizeof(double)))==NULL)

  {
  printf("Erro ao alocar memória");
  return;
  }
  GL(x,w);
}

Agora, implementamos a primeira versão do código que resolve o problema:

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#define M 1000
#define N 2000

const double L=1;
const double B0=1;
const double BL=1;
const double sigma=1;
const double lambda=1.0;
double *mu,*w,*y;

double GL()
{
  size_t i,flag;
  gsl_integration_glfixed_table * table;
  table = gsl_integration_glfixed_table_alloc (2*M);
  for (i=M; i<2*M;i++)
  {
    flag = gsl_integration_glfixed_point (-1, 1, i, &mu[i-M], &w[i-M], table);
    //printf("flag=%zu, %zu %f %f\n",flag, i, mu[i], w[i]);
  }
  gsl_integration_glfixed_table_free(table);
}
double Q(double x)
{
return exp(-x);
}
void main(void)
{

  int i,k,n;
  double **I,*II;
  if (((mu=malloc(sizeof(double)*M))==NULL)||((w=malloc(sizeof(double)*M))==NULL))
  {
  printf("Erro ao alocar memória");
  return;
  }
  GL();
  if ((y=malloc(sizeof(double)*(N+1)))==NULL)
  {
  printf("Erro ao alocar memória");
  return;
  }
  if ((II=malloc(sizeof(double)*(N+1)))==NULL)
  {
  printf("Erro ao alocar memória");
  return;
  }
  if ((I = calloc(N+1,sizeof(double)))==NULL)
  {
  printf("Erro ao alocar memória");
  return;
  }
  for (i=0;i<N+1;i++) if ((I[i] = calloc(2*M,sizeof(double)))==NULL)
  {
  printf("Erro ao alocar memória");
  return;
  }
//malha espacial
  double h = L/N;
  for (i=0;i<N+1;i++) y[i] = i*h;
//Contorno
  for (i=0;i<M;i++)
  {
    I[0][i] = B0;

    I[N][M+i] =BL;
  }
  //iteracao
  double Exp,S=0;
  for (n=0;n<200;n++)
  {
    //Fluxo Escalar
    for(k=0;k<N+1;k++)
    {
      II[k]=0;
      for(i=0;i<M;i++)
      {
        II[k] += w[i]*(I[k][i]+I[k][M+i]);
      }
      if (k%200==0) printf("II(%d)=%f\n",k,II[k]);
    }

    //Recorrência
    for (i=0;i<M;i++)
    {
      for(k=0;k<N;k++)
      {
        S = 1/2.0*(sigma/2.0*(II[k]+II[k+1])+Q(y[k])+Q(y[k+1]));
        Exp = exp(lambda/mu[i]*(y[k]-y[k+1]));

        I[k+1][i] = S/lambda*(1-Exp)+I[k][i]*Exp;
      }

      for(k=N-1;k>-1;k--)
      {
        S = 1/2.0*(sigma/2.0*(II[k]+II[k+1])+Q(y[k])+Q(y[k+1]));
        Exp = exp(lambda/mu[i]*(y[k]-y[k+1]));
        I[k][M+i] = S/lambda*(1-Exp)+I[k+1][M+i]*Exp;

      }

    }
  }
}

Para validar usamos resultados de artigos científicos, que fornecem dados tais como dessa tabela:

Tabela 13.1: Valores calculados para

\int_{- 1}^{1} I (y, μ) d μ

quando

Q (y) = e^{- y}

λ^{'} = σ^{'} = 1

B_{0} = B_{L} = 1

, com

L = 1


y	0.0	0.2	0.4	0.6	0.8	1.0

$G F D_{800}$	3.514742	4.193467	4.307001	4.162773	3.820960	3.196349

DD & Gauss-Legendre	3.514748	4.193470	4.307004	4.162776	3.820963	3.196353

DD & Clenshaw-Curtis	3.514748	4.193471	4.307005	4.162777	3.820963	3.196354

Naturalmente, o código pode ser melhorado. Por exemplo, podemos estabelecer um critério de convergência para o processo iterativo.

Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 19/8/2020 às 17:33:16.

Recursos

Projeto

IME - UFRGS
- Página do IME
- Página da UFRGS

UFRGS - IME - Recursos Educacionais Abertos de Matemática. Contato: reamat@ufrgs.br.

Informe erros ou edite você mesmo!