Problemas

Computação Científica - Um Livro Colaborativo

4.3 Problemas

Exemplo 4.3.1. Resolva numericamente o problema de valor de contorno $\begin{array}{rcl} - u_{x x} & = & 100 {(x - 1)}^{2}, 0 < x < 1, \\ u (0) & = & 0, \\ u (1) & = & 0 . \end{array}$

usando a fórmula de diferenças finitas central de ordem 2 para discretizar a derivada. Compute o erro absoluto médio definido por

E : = \frac{1}{N} \sum_{i = 1}^{N} |u (x_{i}) - u_{i}|,

onde $x_{i}$ é o $i$ -ésimo ponto da malha, $i = 1, 2, \dots, N$ e $N$ é o número de pontos da mesma, $u_{i}$ é a solução aproximada e $u (x_{i})$ é a solução analítica, dada por

u (x) = - \frac{25 {(x - 1)}^{4}}{3} - \frac{25 x}{3} + \frac{25}{3} .

Solução: Primeiro, vamos montar a discretização do problema que está definido no domínio $[0, 1]$ . A malha é dada por:

x_{i} = (i - 1) h, i = 1, 2, \dots, N,

com $h = 1 ∕ (N - 1)$ .

Aplicando na equação diferencial o método de diferenças finitas central de ordem 2 para aproximar a derivada segunda, obtemos as seguintes $N - 2$ equações:

- \frac{u_{i - 1} - 2 u_{i} + u_{i + 1}}{h^{2}} = 100 {(x_{i} - 1)}^{2}, i = 2, \dots, N - 1 .

As equações de contorno levam as condições:

u_{1} = u_{N} = 0 .

Ou seja, obtemos o seguinte sistema linear $N \times N$ : $\begin{array}{rcl} u_{1} & = & 0, & (4.2) \\ - \frac{1}{h^{2}} (u_{i - 1} - 2 u_{i} + u_{i + 1}) & = & 100 {(x_{i} - 1)}^{2}, i = 2, \dots, N - 1, & (4.3) \\ u_{N} & = & 0 . & (4.4) \end{array}$

Observamos que este é um sistema linear $N \times N$ , o qual pode ser escrito na forma matricial $A \underset{}{u} = b$ , cujos matriz de coeficientes é

A = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & \dots & 0 \\ 1 & - 2 & 1 & 0 & 0 & \dots & 0 \\ 0 & 1 & - 2 & 1 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & \dots & 1 \end{matrix}],

o vetor das incógnitas e o vetor dos termos constantes são

\underset{}{u} = [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \\ u_{N} \end{matrix}] e b = [\begin{matrix} 0 \\ - 100 h^{2} {(x_{2} - 1)}^{2} \\ - 100 h^{2} {(x_{3} - 1)}^{2} \\ ⋮ \\ 0 \end{matrix}] .

Vamos aplicar o método de Jacobi para resolver esse problema. Partimos do exemplo 4.2.2 e introduzimos a matriz oriunda da discretização:

#include <stdio.h>
#include <math.h>

#define N 101 /* dimensão do sistema*/

double norma_2(double x[N])
{
  int i;
  double norma2=x[0]*x[0]/(2*(N-1));
  for (i=1;i<N-1;i++)
  {
    norma2+=x[i]*x[i]/(N-1);
  }
  norma2+=x[N-1]*x[N-1]/(2*(N-1));
  return sqrt(norma2);
}
void Iteracao(double u_antigo[N],double u_atual[N],double matriz[N][N], double vetor[N])
{
  double aux;
  int i,j;
  for (i=0;i<N;i++)
  {
    aux=0;
    for (j=0;j<i;j++)
    aux+=matriz[i][j]*u_antigo[j];
    for (j=i+1;j<N;j++)
    aux+=matriz[i][j]*u_antigo[j];
    u_atual[i]=(vetor[i]-aux)/matriz[i][i];
  }
}
void Jacobi(double tolerancia,double u_atual[N],double vetor[N],double matriz[N][N],double malha[N])
{

  double u_antigo[N],dif[N];
  int i,controle=3;

  //chute inicial
  for (i=0;i<N;i++)u_antigo[i]=0;

  //iteração
  while (controle)
  {
    Iteracao(u_antigo,u_atual,matriz,vetor);
    for (i=0;i<N;i++) dif[i]=u_atual[i]-u_antigo[i];
    if (norma_2(dif)<tolerancia) controle--;
    else controle=3;
    for (i=0;i<N;i++) u_antigo[i]=u_atual[i];
  }
  return;
}

double solucao_analitica(double x)
{
return -25.*(x-1.)*(x-1.)*(x-1.)*(x-1.)/3-25.*x/3+25./3.;
}

main (void)
{
  int i,j;
  double tolerancia=1e-5,solucao[N];
  double matriz[N][N];
  double vetor[N],x[N];
  //malha
  double h=1./(N-1);
  for (i=0;i<N;i++) x[i]=i*h;

  //matriz e vetor
  for (i=0;i<N;i++) for (j=0;j<N;j++) matriz[i][j]=0;
  matriz[0][0]=1;
  vetor[0]=0;
  for (i=1;i<N-1;i++)
  {
    matriz[i][i]=-2;
    matriz[i][i+1]=1;
    matriz[i][i-1]=1;
    vetor[i]=-100*h*h*(x[i]-1)*(x[i]-1);
  }
  matriz[N-1][N-1]=1;
  vetor[N-1]=0;

  Jacobi(tolerancia,solucao,vetor,matriz,x);
  //erro médio
  double erro=0;
  for (i=0;i<N;i++) erro+=fabs(solucao[i]-solucao_analitica(x[i]));
  erro/=N;
  printf("erro médio = %f\n",erro);
  //for (i=0;i<N;i++)     printf("u_%d=%f, u(x[%d]) = %f\n",i,solucao[i],i,
  solucao_analitica(x[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!