Problemas

Computação Científica - Um Livro Colaborativo

6.3 Problemas

Você deve ter observado que todos os exercícios na aula 8 resultam em um sistema tridiagonal. O algorítmo de Thomas é um método eficiente para resolvê-los.

Exemplo 6.3.1. Implemente o método de Thomas para resolver um sistema $A x = b$ de dimensão $N$ .

O algoritmo TDMA para resolver o sistema $\begin{array}{rcl} [\begin{matrix} b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1} \\ d_{2} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}] . \end{array}$

é dado por $\begin{array}{rcl} c_{i}^{'} & = & \{\begin{matrix} \begin{matrix} \frac{c_{i}}{b_{i}}, & i = 1 \\ \frac{c_{i}}{b_{i} - a_{i} c_{i - 1}^{'}}, & i = 2, 3, \dots, n - 1 \end{matrix} \end{matrix}, \\ d_{i}^{'} & = & \{\begin{matrix} \begin{matrix} \frac{d_{i}}{b_{i}}, & i = 1 \\ \frac{d_{i} - a_{i} d_{i - 1}^{'}}{b_{i} - a_{i} c_{i - 1}^{'}}, & i = 2, 3, \dots, n . \end{matrix} \end{matrix} \end{array}$

e $\begin{array}{rcl} x_{n} & = & d_{n}^{'} \\ x_{i} & = & d_{i}^{'} - c_{i}^{'} x_{i + 1}, i = n - 1, n - 2, \dots, 1 . \end{array}$

#include <stdio.h>
#include <math.h>
#define N 3

//Método de Thomas para resolver Ax=y
//x entra o vetor y e saí a solução x (N posições, i=0,....,N-1)
//a entra a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b entra a diagonal  (N posições, i=0,....,N-1)
//c entra a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void Thomas(double x[], const double a[], const double b[], double c[])
{
  int i;

  /*calculo de c[0]’ e d[0]’*/
  c[0] = c[0] / b[0];
  x[0] = x[0] / b[0];

  /* laço para calcular c’ e d’ */
  for (i = 1; i < N; i++)
  {
    double aux = 1.0/ (b[i] - a[i] * c[i - 1]);
    c[i] = c[i] * aux;
    x[i] = (x[i] - a[i] * x[i - 1]) * aux;
  }

  /* Calculando a solução */
  for (i = N - 1; i >= 0; i--)
    x[i] = x[i] - c[i] * x[i + 1];
}

main (void)
{
  double x[N],a[N],b[N],c[N];
  a[0]=0;c[N-1]=0;

  x[0]=4;x[1]=9;x[2]=11;
  a[1]=1;a[2]=1;
  b[0]=3;b[1]=3;b[2]=3;
  c[0]=1;c[1]=1;
  Thomas(x,a,b,c);
  printf("%f, %f, %f\n",x[0],x[1],x[2]);
}

Exemplo 6.3.2. Use o algoritmo TDMA para resolver o PVC $\{\begin{matrix} - u_{x x} + u_{x} = 200 e^{- {(x - 1)}^{2}}, 0 < x < 1 . \\ 15 u (0) + u^{'} (0) = 500 \\ 10 u (1) + u^{'} (1) = 1 \end{matrix}$

Uma versão discreta oriunda do método de diferenças finitas é $\begin{array}{l} 15 u_{1} + \frac{u_{2} - u_{1}}{h} = 500 \\ \frac{- u_{i + 1} + 2 u_{i} - u_{i - 1}}{h^{2}} + \frac{u_{i + 1} - u_{i - 1}}{2 h} = 200 e^{- {(x_{i} - 1)}^{2}}, i = 2, 3, . . ., N - 1 \\ 10 u_{N} + \frac{u_{N} - u_{N - 1}}{h} = 1 \end{array}$

onde $x_{i} = h (i - 1)$ , $i = 1, 2, . . ., N$ , $h = \frac{1}{N - 1}$ . Isto é, $\begin{array}{l} u_{1} (15 - \frac{1}{h}) + u_{2} (\frac{1}{h}) = 500 \\ u_{i - 1} (- \frac{1}{h^{2}} - \frac{1}{2 h}) + u_{i} (\frac{2}{h^{2}}) + u_{i + 1} (- \frac{1}{h^{2}} + \frac{1}{2 h}) = 200 e^{- {(x_{i} - 1)}^{2}}, i = 2, 3, . . ., N - 1 \\ u_{N - 1} (- \frac{1}{h}) + u_{N} (10 + \frac{1}{h}) = 1 \end{array}$

ou ainda

[\begin{matrix} 15 - \frac{1}{h} & \frac{1}{h} & 0 & 0 & 0 & \dots & 0 \\ - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} & 0 & 0 & \dots & 0 \\ 0 & - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 & 0 & - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} \\ 0 & \dots & 0 & 0 & 0 & - \frac{1}{h} & 10 + \frac{1}{h} \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \\ u_{N - 1} \\ u_{N} \end{matrix}] = [\begin{matrix} 500 \\ 200 e^{- {(x_{2} - 1)}^{2}} \\ 200 e^{- {(x_{3} - 1)}^{2}} \\ ⋮ \\ 200 e^{- {(x_{N - 1} - 1)}^{2}} \\ 1 \end{matrix}]

#include <stdio.h>
#include <math.h>
#define N 11

//Método de Thomas para resolver Ax=y
//x entra o vetor y e saí a solução x (N posições, i=0,....,N-1)
//a entra a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b entra a diagonal  (N posições, i=0,....,N-1)
//c entra a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void Thomas(double x[], const double a[], const double b[], double c[])
{
  int i;

  /*calculo de c[0]’ e d[0]’*/
  c[0] = c[0] / b[0];
  x[0] = x[0] / b[0];

  /* laço para calcular c’ e d’ */
  for (i = 1; i < N; i++)
  {
    double aux = 1.0/ (b[i] - a[i] * c[i - 1]);
    c[i] = c[i] * aux;
    x[i] = (x[i] - a[i] * x[i - 1]) * aux;
  }

  /* Calculando a solução */
  for (i = N - 1; i >= 0; i--)
    x[i] = x[i] - c[i] * x[i + 1];
}

main (void)
{
  double x[N],a[N],b[N],c[N];
  a[0]=0;c[N-1]=0;

  //malha
  int i;
  double h=1./(N-1);
  double p[N];
  for (i=0;i<N;i++) p[i]=i*h;

  //sistema
  b[0]=15-1/(h);
  c[0]=1/(h);
  x[0]=500;
  for (i=1;i<N-1;i++)
  {
  b[i]=2/(h*h);
  a[i]=-1/(h*h)-1/(2*h);
  c[i]=-1/(h*h)+1/(2*h);
  x[i]=200*exp(-(p[i]-1)*(p[i]-1));
  }
  b[N-1]=10+1/h;
  a[N-1]=-1/h;
  x[N-1]=1;
  Thomas(x,a,b,c);
  for (i=0;i<N;i++) printf("%f\n",x[i]);
}

Com o objetivo de resolver alguns problemas de contorno não lineares, vamos implementar o método de Newton para sistema.

Exemplo 6.3.3. Implemente um código para calcular raízes do sistema não linear $\begin{array}{rcl} x_{1}^{2} = cos (x_{1} x_{2}) + 1 \\ sen (x_{2}) = 2 cos (x_{1}) . \end{array}$

Calcule uma aproximação para a raíz que fica próxima ao ponto $x_{1} = 1, 5$ e $x_{2} = 0, 5$ .

O método iterativo de Newton-Raphson para encontrar as raízes de $F (x) = 0$ , $F (x) = [\begin{matrix} f_{1} (x_{1}, x_{2}, \dots, x_{n}) \\ f_{2} (x_{1}, x_{2}, \dots, x_{n}) \\ ⋮ \\ f_{n} (x_{1}, x_{2}, \dots, x_{n}) \end{matrix}],$

é dado por

\{\begin{matrix} x^{(0)} & = & dado inicial \\ x^{(k + 1)} & = & x^{(k)} - J_{F}^{- 1} (x^{(k)}) F (x^{(k)}), k \geq 0 \end{matrix}

onde a matriz Jacobina $J_{F}$ é $J_{F} = \frac{\partial (f_{1}, f_{2}, \dots, f_{n})}{\partial (x_{1}, x_{2}, \dots, x_{n})} = [\begin{matrix} \frac{\partial f_{1}}{\partial x_{1}} & \frac{\partial f_{1}}{\partial x_{2}} & \dots & \frac{\partial f_{1}}{\partial x_{n}} \\ \frac{\partial f_{2}}{\partial x_{1}} & \frac{\partial f_{2}}{\partial x_{2}} & \dots & \frac{\partial f_{2}}{\partial x_{n}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{n}}{\partial x_{1}} & \frac{\partial f_{n}}{\partial x_{2}} & \dots & \frac{\partial f_{n}}{\partial x_{n}} \end{matrix}] .$

Especificamente no problema do exemplo 6.3.3, temos $\begin{array}{rcl} F (x) = [\begin{matrix} x_{1}^{2} - cos (x_{1} x_{2}) - 1 \\ sen (x_{2}) - 2 cos (x_{1}) \end{matrix}] . \end{array}$

e $\begin{array}{rcl} J_{F} (x) = [\begin{matrix} 2 x_{1} + x_{2} sen (x_{1} x_{2}) & x_{1} sen (x_{1} x_{2}) \\ 2 sen (x_{1}) & cos (x_{2}) \end{matrix}] . \end{array}$

Observamos que na prática não invertemos a matriz Jacobina, mas calculamos a solução do sistema $J_{F} δ = F$ :

\{\begin{matrix} x^{(0)} & = & [\begin{matrix} 1, 5 \\ 0, 5 \end{matrix}] \\ J_{F} (x^{(k)}) δ^{(k)} & = & F (x^{(k)}) (Calcula-se a solução δ^{(k)} do sistema) \\ x^{(k + 1)} & = & x^{(k)} - δ^{(k)}, k \geq 0 \end{matrix}

#include <stdio.h>
#include <math.h>
#define N 2

//Método de Thomas para resolver Ax=y
//x entra o vetor y e saí a solução x (N posições, i=0,....,N-1)
//a entra a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b entra a diagonal  (N posições, i=0,....,N-1)
//c entra a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void Thomas(double x[], const double a[], const double b[], double c[])
{
  int i;

  /*calculo de c[0]’ e d[0]’*/
  c[0] = c[0] / b[0];
  x[0] = x[0] / b[0];

  /* laço para calcular c’ e d’ */
  for (i = 1; i < N; i++)
  {
    double aux = 1.0/ (b[i] - a[i] * c[i - 1]);
    c[i] = c[i] * aux;
    x[i] = (x[i] - a[i] * x[i - 1]) * aux;
  }

  /* Calculando a solução */
  for (i = N - 1; i >= 0; i--)
    x[i] = x[i] - c[i] * x[i + 1];
}
//Função
//x é o ponto de entrada
//f é o vetor com o resultado de f(x)
void F(double *x,double *f)
{
  f[0]=x[0]*x[0]-cos(x[0]*x[1])-1.0;
  f[1]=sin(x[1])-2*cos(x[0]);
}
//Jacobiana
//x é o ponto de entrada
//a sai a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b sai a diagonal  (N posições, i=0,....,N-1)
//c sai a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void J(double *x,double *a,double *b,double *c)
{
  b[0]=2*x[0]+x[1]*sin(x[0]*x[1]);
  a[0]=0;
  c[0]=x[0]*sin(x[0]*x[1]);
  a[1]=2*sin(x[0]);
  c[1]=0;
  b[1]=cos(x[1]);
}
//norma da diferença entre dois vetores
double norma_dif(double x[N],double y[N])
{
  int i;
  double norma2=0;
  for (i=0;i<N;i++)
  {
    norma2+=(x[i]-y[i])*(x[i]-y[i]);
  }
  return sqrt(norma2);
}
main (void)
{
  int i;
  double x0[N],x[N],d[N],a[N],b[N],c[N],tol=1e-10,*px0=x0,*px=x,*paux=NULL;
  x0[0]=1.5;
  x0[1]=0.5;
  do
  {
    F(px0,d);
    J(px0,a,b,c);
    Thomas(d,a,b,c);
    for (i=0;i<N;i++) *(px+i)=*(px0+i)-d[i];
    for (i=0;i<N;i++) printf("%.12f  ",*(px+i));
    printf("\n");
    paux=px0;
    px0=px;
    px=paux;
  }while(norma_dif(px0,px)>tol);
}

Exemplo 6.3.4. Use o método de diferenças finitas e o método de Newton para resolver o problema de valor de contono não-linear $\{\begin{matrix} - u_{x x} + u_{x} = e^{- u^{2}}, 0 < x < 1 . \\ u^{'} (0) + 3 u (0) = 10 \\ u (1) = 5 \end{matrix}$

Uma versão discreta oriunda do método de diferenças finitas é $\begin{array}{l} 3 u_{1} + \frac{u_{2} - u_{1}}{h} = 10 \\ \frac{- u_{i + 1} + 2 u_{i} - u_{i - 1}}{h^{2}} + \frac{u_{i + 1} - u_{i - 1}}{2 h} = e^{- u_{i}^{2}}, i = 2, 3, . . ., N - 1 \\ u_{N} = 5 \end{array}$

onde $x_{i} = h (i - 1)$ , $i = 1, 2, . . ., N$ , $h = \frac{1}{N - 1}$ . Podemos colocar na forma $F (u) = 0$ , onde $u = {[u_{1} u_{2} \dots u_{N}]}^{T}$ e

F (u) = \{\begin{matrix} 3 u_{1} + \frac{u_{2} - u_{1}}{h} - 10 \\ \frac{- u_{i + 1} + 2 u_{i} - u_{i - 1}}{h^{2}} + \frac{u_{i + 1} - u_{i - 1}}{2 h} - e^{- u_{i}^{2}}, i = 2, 3, . . ., N - 1 \\ u_{N} - 5 \end{matrix}

A jacobina desse sistema não linear é

J_{F} (u) = [\begin{matrix} 3 - \frac{1}{h} & \frac{1}{h} & 0 & 0 & 0 & \dots & 0 \\ - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} + 2 u_{i} e^{- u_{i}^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} & 0 & 0 & \dots & 0 \\ 0 & - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} + 2 u_{i} e^{- u_{i}^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 & 0 & - \frac{1}{h^{2}} - \frac{1}{2 h} & \frac{2}{h^{2}} + 2 u_{i} e^{- u_{i}^{2}} & - \frac{1}{h^{2}} + \frac{1}{2 h} \\ 0 & \dots & 0 & 0 & 0 & 0 & 1 \end{matrix}]

Vamos implementar um código usando o método de Thomas para resolver os sistemas tridiagonais.

#include <stdio.h>
#include <math.h>
#define N 11
#define h (1./(N-1))

//Método de Thomas para resolver Ax=y
//x entra o vetor y e saí a solução x (N posições, i=0,....,N-1)
//a entra a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b entra a diagonal  (N posições, i=0,....,N-1)
//c entra a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void Thomas(double x[], const double a[], const double b[], double c[])
{
  int i;

  /*calculo de c[0]’ e d[0]’*/
  c[0] = c[0] / b[0];
  x[0] = x[0] / b[0];

  /* laço para calcular c’ e d’ */
  for (i = 1; i < N; i++)
  {
    double aux = 1.0/ (b[i] - a[i] * c[i - 1]);
    c[i] = c[i] * aux;
    x[i] = (x[i] - a[i] * x[i - 1]) * aux;
  }

  /* Calculando a solução */
  for (i = N - 1; i >= 0; i--)
    x[i] = x[i] - c[i] * x[i + 1];
}
//Função
//x é o ponto de entrada
//f é o vetor com o resultado de f(x)
void F(double *x,double *f)
{
  int i;
  f[0]=3*x[0]+(x[1]-x[0])/h-10;
  for (i=1;i<N-1;i++) f[i]=(-x[i+1]+2*x[i]-x[i-1])/(h*h)+(x[i+1]-x[i-1])/(2*h)-exp(-x[i]*x[i]);
  f[N-1]=x[N-1]-5;
}
//Jacobiana
//x é o ponto de entrada
//a sai a subdiagonal à esquerda de A (N-1 posições, i=1,...,N-1)
//b sai a diagonal  (N posições, i=0,....,N-1)
//c sai a subdiagonal à direita de A (N-1 posições, i=0,...,N-2)
void J(double *x,double *a,double *b,double *c)
{
  int i;
  a[0]=0;
  b[0]=3.-1./h;
  c[0]=1./h;
  for (i=1;i<N-1;i++)
  {
    a[i]=-1./(h*h)-1./(2*h);
    b[i]=2./(h*h)+2*x[i]*exp(-x[i]*x[i]);
    c[i]=-1./(h*h)+1./(2*h);
  }
  a[N-1]=0;
  b[N-1]=1.;
  c[N-1]=0;
}
//norma da diferença entre dois vetores
double norma_max(double x[N])
{
  int i;
  double norm=fabs(x[0]);
  for (i=0;i<N;i++) if(norm<fabs(x[i])) norm=fabs(x[i]);
  return norm;
}
void main (void)
{
  int i;
  double x0[N],x[N],d[N],a[N],b[N],c[N],tol=1e-10,*px0=x0,*px=x,*paux=NULL;
  for (i=0;i<N;i++) x0[i]=0;
  do
  {
    F(px0,d);
    J(px0,a,b,c);
    Thomas(d,a,b,c);
    for (i=0;i<N;i++) *(px+i)=*(px0+i)-d[i];
    for (i=0;i<N;i++) printf("%.12f  ",*(px+i));
    printf("\n");
    paux=px0;
    px0=px;
    px=paux;
  }while(norma_max(d)>tol);
  for (i=0;i<N;i++) printf("%.12f  ",*(px+i));
}

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!