Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

| | | |

Matemática Numérica Paralela

2 Multiprocessamento (MP)2.5 Aplicação: Equação de Poisson 2.7 Aplicação: Equação de Fisher

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

2.6 Aplicação: Equação do Calor

Em construção

Vamos considerar a equação do Calor com dadas condições inicial e de contorno

	$\displaystyle u_{t}=\alpha u_{xx}+f,\leavevmode\nobreak\ (t,x)\in(0,t_{f})% \times(a,b),$		(2.4a)
	$\displaystyle u(0,x)=u^{0}(x),\leavevmode\nobreak\ x\in[a,b],$		(2.4b)
	$\displaystyle u(t,a)=u(t,b)=0,$		(2.4c)

com dado $\alpha>0$ e fonte $f$ .

2.6.1 Formulação Explícita

Assumimos uma discretização no tempo com $n_{t}$ passos $t^{(k)}=kh_{t}$ , de tamanho $h_{t}=t_{f}/n_{t}$ , $k=0,1,\dotsc,n_{t}$ . Denotando $u^{(k)}(x)\approx u\left(t^{(k)},x\right)$ , temos o esquema de Euler explícito

$u^{(0)}=u^{0}(x),$		(2.5a)
	$\displaystyle u^{(k+1)}=u^{(k)}+h_{t}\left[\alpha u^{(k)}_{xx}+f^{(k)}(x)% \right],$	(2.5b)
	$\displaystyle u^{(k+1)}(a)=u^{(k+1)}(t,b)=0.$	(2.5b)

Agora, assumimos uma malha espacial com $n_{x}+1$ nodos $x_{i}=a+ih_{x}$ , de tamanho $h_{x}=(b-a)/n_{x}$ , $i=0,1,\dotsc,n_{x}$ , e denotando $u^{(k)}_{i}\approx u\left(t^{(k)},x_{i}\right)$ e usando a fórmula de diferenças finitas central, obtemos o esquema iterativo

	$\displaystyle u^{(0)}_{i}=u^{0}(x_{i}),\leavevmode\nobreak\ \forall i,$		(2.6a)
	$\displaystyle u^{(k+1)}_{i}=u^{(k)}_{i}+ru^{(k)}_{i-1}-2ru^{(k)}_{i}+ru^{(k)}_% {i+1}+h_{t}f^{(k)}_{i},$		(2.6b)
	$\displaystyle u^{(k+1)}_{0}=u^{(k+1)}_{n_{x}}=0,$		(2.6c)

onde $r:=\alpha h_{t}/h_{x}^{2}$ .

Observação 2.6.1.

O esquema de Euler explícito (2.6.1) requer a seguinte condição de estabilidade

\frac{h_{t}}{h_{x}^{2}}\leq\frac{1}{2}.

(2.7)

Exemplo 2.6.1.

Consideramos o seguinte problema de calor

	$\displaystyle u_{t}=u_{xx}+(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ (t,x)\in(0,1)^{2},$		(2.8a)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(2.8b)
	$\displaystyle u(t,0)=u(t,1)=0.$		(2.8c)

Código 3: ex_mlp_calor_explicito.cc

⬇

1 #include <stdio.h>

2 #include <stdlib.h>

4 #include <omp.h>

6 #include <math.h>

7 #include <gsl/gsl_vector.h>

9 // fonte

10 double f(double t, double x) {

11 return (M_PI*M_PI-1.) *

12 exp(-t) * sin(M_PI * x);

13 }

15 int main()

16 {

17 // parâmetros no tempo

18 double tf = 1.;

19 int nt = 1000;

20 double ht = tf / nt;

22 // parâmetros no espaço

23 int nx = 10;

24 double hx = 1./nx;

26 // malha

27 gsl_vector *x = gsl_vector_alloc(nx+1);

28 #pragma omp parallel for

29 for (int i = 0; i <= nx; i++) {

30 gsl_vector_set(x, i, i*hx);

31 }

33 gsl_vector *u0 = gsl_vector_alloc(nx+1);

34 gsl_vector *u = gsl_vector_calloc(nx+1);

36 // inicialização

37 double t = 0.;

38 #pragma omp parallel for

39 for (int i = 0; i <= nx; i++) {

40 gsl_vector_set(u0, i,

41 sin(M_PI *

42 gsl_vector_get(x, i)));

43 }

45 // armazena malha

46 FILE *fx = fopen("results/x.dat", "wb");

47 gsl_vector_fwrite(fx, x);

48 fclose(fx);

50 // armazena solução

51 char str[21];

52 sprintf(str, "results/u_%06d.dat", 0);

53 printf("t = %g, [%s]\n", t, str);

55 FILE *fu = fopen(str, "wb");

56 gsl_vector_fwrite(fu, u0);

57 fclose(fu);

59 #pragma omp parallel

60 {

61 // iteração no tempo

62 for (int k=0; k < nt; k++) {

63 #pragma omp single

64 t += ht;

66 #pragma omp for

67 for (int i=1; i < nx; i++) {

68 gsl_vector_set(u, i,

69 gsl_vector_get(u0, i) +

70 ht/(hx*hx) * (gsl_vector_get(u0, i+1) - 2*gsl_vector_get(u0, i) + gsl_vector_get(u0, i-1)) + ht *

71 f(t-ht, gsl_vector_get(x, i)));

72 }

74 #pragma omp single

75 {

76 sprintf(str, "results/u_%06d.dat", k+1);

77 printf("t = %g, [%s]\n", t, str);

78 fu = fopen(str, "wb");

79 gsl_vector_fwrite(fu, u0);

80 fclose(fu);

82 // u0 = u

83 gsl_vector_memcpy(u0, u);

84 }

86 }

87 }

89 gsl_vector_free(x);

90 gsl_vector_free(u);

91 gsl_vector_free(u0);

92 return 0;

93 }

2.6.2 Formulação Implícita

A condição de estabilidade (2.7) pode ser contornada aplicando-se o esquema de Euler implícito

	$\displaystyle u^{(0)}_{i}=u^{0}(x_{i}),\leavevmode\nobreak\ \forall i,$		(2.9a)
	$\displaystyle-ru^{(k+1)}_{i-1}+(2r+1)u^{(k+1)}_{i}-ru^{(k+1)}_{i+1}=u^{(k)}_{i% }+h_{t}f^{(k+1)}_{i},$		(2.9b)
	$\displaystyle u^{(k+1)}_{0}=u^{(k+1)}_{n_{x}}=0,$		(2.9c)

onde $r:=\alpha h_{t}/h_{x}^{2}$ .

Exemplo 2.6.2.

Implementamos o esquema de Euler implícito (2.6.2) para o problema de calor (2.6.1) do exemplo anterior. No código abaixo, o sistema linear em cada iteração é computado com o método de Jacobi.

Código 4: ex_mlp_calor_implicito.cc

⬇

1#include <stdio.h>

2#include <stdlib.h>

4#include <omp.h>

6#include <math.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_blas.h>

10// fonte

11double f(double t, double x) {

12 return (M_PI*M_PI-1.) *

13 exp(-t) * sin(M_PI * x);

14}

16int main()

17{

18 // parâmetros no tempo

19 double tf = 1.;

20 int nt = 100;

21 double ht = tf / nt;

23 // parâmetros no espaço

24 int nx = 10;

25 double hx = 1./nx;

27 // parâmetros Jacobi

28 int maxit = 1000;

29 double tol = 1.49e-8;

31 // auxiliar

32 double r = ht/(hx*hx);

34 // malha

35 gsl_vector *x = gsl_vector_alloc(nx+1);

36 #pragma omp parallel for

37 for (int i = 0; i <= nx; i++) {

38 gsl_vector_set(x, i, i*hx);

39 }

41 gsl_vector *u0 = gsl_vector_alloc(nx+1);

42 gsl_vector *u = gsl_vector_calloc(nx+1);

44 // auxiliares

45 double norm;

46 gsl_vector *u00 = gsl_vector_alloc(nx+1);

47 gsl_vector *diff = gsl_vector_calloc(nx+1);

49 // inicialização

50 double t = 0.;

51 #pragma omp parallel for

52 for (int i = 0; i <= nx; i++) {

53 gsl_vector_set(u0, i,

54 sin(M_PI *

55 gsl_vector_get(x, i)));

56 }

58 // u00 = u0

59 gsl_vector_memcpy(u00, u0);

61 // armazena malha

62 FILE *fx = fopen("results/x.dat", "wb");

63 gsl_vector_fwrite(fx, x);

64 fclose(fx);

66 // armazena solução

67 char str[21];

68 sprintf(str, "results/u_%06d.dat", 0);

69 printf("t = %g, [%s]\n", t, str);

71 FILE *fu = fopen(str, "wb");

72 gsl_vector_fwrite(fu, u0);

73 fclose(fu);

75 #pragma omp parallel

76 {

77 // iteração no tempo

78 for (int k=0; k < nt; k++) {

80 #pragma omp single

81 {

82 t += ht;

83 printf("%d: t = %g\n", k+1, t);

84 }

86 // iteração de Jacobi

87 for (int it=0; it < maxit; it++) {

89 #pragma omp for

90 for (int i=1; i < nx; i++) {

91 gsl_vector_set(u, i, 1./(2.*r + 1.) *

92 (r * (gsl_vector_get(u00, i-1) +

93 gsl_vector_get(u00, i+1)) +

94 gsl_vector_get(u0, i) +

95 + ht * f(t, gsl_vector_get(x, i))));

97 gsl_vector_set(diff, i, gsl_vector_get(u, i) -

98 gsl_vector_get(u00, i));

99 }

100

101 #pragma omp single

102 {

103 norm = gsl_blas_dnrm2(diff);

104 printf("\t%d: norm=%.2e\n", it, norm);

105

106 gsl_vector_memcpy(u00, u);

107 }

108

109 if (norm < tol) {

110 break;

111 }

112

113 }

114

115 #pragma omp single

116 {

117 sprintf(str, "results/u_%06d.dat", k+1);

118 printf("t = %g, [%s]\n", t, str);

119 fu = fopen(str, "wb");

120 gsl_vector_fwrite(fu, u0);

121 fclose(fu);

122

123 // u0 = u

124 gsl_vector_memcpy(u0, u);

125 }

126

127 }

128

129 }

130

131 gsl_vector_free(x);

132 gsl_vector_free(u);

133 gsl_vector_free(u0);

134 gsl_vector_free(u00);

135 gsl_vector_free(diff);

136

137 return 0;

138}

2.6.3 Exercícios

Em construção

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Matemática Numérica Paralela

2 Multiprocessamento (MP)2.5 Aplicação: Equação de Poisson 2.7 Aplicação: Equação de Fisher

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

2.6 Aplicação: Equação do Calor

Em construção

Vamos considerar a equação do Calor com dadas condições inicial e de contorno

	$\displaystyle u_{t}=\alpha u_{xx}+f,\leavevmode\nobreak\ (t,x)\in(0,t_{f})% \times(a,b),$		(2.4a)
	$\displaystyle u(0,x)=u^{0}(x),\leavevmode\nobreak\ x\in[a,b],$		(2.4b)
	$\displaystyle u(t,a)=u(t,b)=0,$		(2.4c)

com dado $\alpha>0$ e fonte $f$ .

2.6.1 Formulação Explícita

$u^{(0)}=u^{0}(x),$		(2.5a)
	$\displaystyle u^{(k+1)}=u^{(k)}+h_{t}\left[\alpha u^{(k)}_{xx}+f^{(k)}(x)% \right],$	(2.5b)
	$\displaystyle u^{(k+1)}(a)=u^{(k+1)}(t,b)=0.$	(2.5b)

	$\displaystyle u^{(0)}_{i}=u^{0}(x_{i}),\leavevmode\nobreak\ \forall i,$		(2.6a)
	$\displaystyle u^{(k+1)}_{i}=u^{(k)}_{i}+ru^{(k)}_{i-1}-2ru^{(k)}_{i}+ru^{(k)}_% {i+1}+h_{t}f^{(k)}_{i},$		(2.6b)
	$\displaystyle u^{(k+1)}_{0}=u^{(k+1)}_{n_{x}}=0,$		(2.6c)

onde $r:=\alpha h_{t}/h_{x}^{2}$ .

Observação 2.6.1.

O esquema de Euler explícito (2.6.1) requer a seguinte condição de estabilidade

\frac{h_{t}}{h_{x}^{2}}\leq\frac{1}{2}.

(2.7)

Exemplo 2.6.1.

Consideramos o seguinte problema de calor

	$\displaystyle u_{t}=u_{xx}+(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ (t,x)\in(0,1)^{2},$		(2.8a)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(2.8b)
	$\displaystyle u(t,0)=u(t,1)=0.$		(2.8c)

Código 3: ex_mlp_calor_explicito.cc

⬇

1 #include <stdio.h>

2 #include <stdlib.h>

4 #include <omp.h>

6 #include <math.h>

7 #include <gsl/gsl_vector.h>

9 // fonte

10 double f(double t, double x) {

11 return (M_PI*M_PI-1.) *

12 exp(-t) * sin(M_PI * x);

13 }

15 int main()

16 {

17 // parâmetros no tempo

18 double tf = 1.;

19 int nt = 1000;

20 double ht = tf / nt;

22 // parâmetros no espaço

23 int nx = 10;

24 double hx = 1./nx;

26 // malha

27 gsl_vector *x = gsl_vector_alloc(nx+1);

28 #pragma omp parallel for

29 for (int i = 0; i <= nx; i++) {

30 gsl_vector_set(x, i, i*hx);

31 }

33 gsl_vector *u0 = gsl_vector_alloc(nx+1);

34 gsl_vector *u = gsl_vector_calloc(nx+1);

36 // inicialização

37 double t = 0.;

38 #pragma omp parallel for

39 for (int i = 0; i <= nx; i++) {

40 gsl_vector_set(u0, i,

41 sin(M_PI *

42 gsl_vector_get(x, i)));

43 }

45 // armazena malha

46 FILE *fx = fopen("results/x.dat", "wb");

47 gsl_vector_fwrite(fx, x);

48 fclose(fx);

50 // armazena solução

51 char str[21];

52 sprintf(str, "results/u_%06d.dat", 0);

53 printf("t = %g, [%s]\n", t, str);

55 FILE *fu = fopen(str, "wb");

56 gsl_vector_fwrite(fu, u0);

57 fclose(fu);

59 #pragma omp parallel

60 {

61 // iteração no tempo

62 for (int k=0; k < nt; k++) {

63 #pragma omp single

64 t += ht;

66 #pragma omp for

67 for (int i=1; i < nx; i++) {

68 gsl_vector_set(u, i,

69 gsl_vector_get(u0, i) +

70 ht/(hx*hx) * (gsl_vector_get(u0, i+1) - 2*gsl_vector_get(u0, i) + gsl_vector_get(u0, i-1)) + ht *

71 f(t-ht, gsl_vector_get(x, i)));

72 }

74 #pragma omp single

75 {

76 sprintf(str, "results/u_%06d.dat", k+1);

77 printf("t = %g, [%s]\n", t, str);

78 fu = fopen(str, "wb");

79 gsl_vector_fwrite(fu, u0);

80 fclose(fu);

82 // u0 = u

83 gsl_vector_memcpy(u0, u);

84 }

86 }

87 }

89 gsl_vector_free(x);

90 gsl_vector_free(u);

91 gsl_vector_free(u0);

92 return 0;

93 }

2.6.2 Formulação Implícita

A condição de estabilidade (2.7) pode ser contornada aplicando-se o esquema de Euler implícito

	$\displaystyle u^{(0)}_{i}=u^{0}(x_{i}),\leavevmode\nobreak\ \forall i,$		(2.9a)
	$\displaystyle-ru^{(k+1)}_{i-1}+(2r+1)u^{(k+1)}_{i}-ru^{(k+1)}_{i+1}=u^{(k)}_{i% }+h_{t}f^{(k+1)}_{i},$		(2.9b)
	$\displaystyle u^{(k+1)}_{0}=u^{(k+1)}_{n_{x}}=0,$		(2.9c)

onde $r:=\alpha h_{t}/h_{x}^{2}$ .

Exemplo 2.6.2.

Implementamos o esquema de Euler implícito (2.6.2) para o problema de calor (2.6.1) do exemplo anterior. No código abaixo, o sistema linear em cada iteração é computado com o método de Jacobi.

Código 4: ex_mlp_calor_implicito.cc

⬇

1#include <stdio.h>

2#include <stdlib.h>

4#include <omp.h>

6#include <math.h>

7#include <gsl/gsl_vector.h>

8#include <gsl/gsl_blas.h>

10// fonte

11double f(double t, double x) {

12 return (M_PI*M_PI-1.) *

13 exp(-t) * sin(M_PI * x);

14}

16int main()

17{

18 // parâmetros no tempo

19 double tf = 1.;

20 int nt = 100;

21 double ht = tf / nt;

23 // parâmetros no espaço

24 int nx = 10;

25 double hx = 1./nx;

27 // parâmetros Jacobi

28 int maxit = 1000;

29 double tol = 1.49e-8;

31 // auxiliar

32 double r = ht/(hx*hx);

34 // malha

35 gsl_vector *x = gsl_vector_alloc(nx+1);

36 #pragma omp parallel for

37 for (int i = 0; i <= nx; i++) {

38 gsl_vector_set(x, i, i*hx);

39 }

41 gsl_vector *u0 = gsl_vector_alloc(nx+1);

42 gsl_vector *u = gsl_vector_calloc(nx+1);

44 // auxiliares

45 double norm;

46 gsl_vector *u00 = gsl_vector_alloc(nx+1);

47 gsl_vector *diff = gsl_vector_calloc(nx+1);

49 // inicialização

50 double t = 0.;

51 #pragma omp parallel for

52 for (int i = 0; i <= nx; i++) {

53 gsl_vector_set(u0, i,

54 sin(M_PI *

55 gsl_vector_get(x, i)));

56 }

58 // u00 = u0

59 gsl_vector_memcpy(u00, u0);

61 // armazena malha

62 FILE *fx = fopen("results/x.dat", "wb");

63 gsl_vector_fwrite(fx, x);

64 fclose(fx);

66 // armazena solução

67 char str[21];

68 sprintf(str, "results/u_%06d.dat", 0);

69 printf("t = %g, [%s]\n", t, str);

71 FILE *fu = fopen(str, "wb");

72 gsl_vector_fwrite(fu, u0);

73 fclose(fu);

75 #pragma omp parallel

76 {

77 // iteração no tempo

78 for (int k=0; k < nt; k++) {

80 #pragma omp single

81 {

82 t += ht;

83 printf("%d: t = %g\n", k+1, t);

84 }

86 // iteração de Jacobi

87 for (int it=0; it < maxit; it++) {

89 #pragma omp for

90 for (int i=1; i < nx; i++) {

91 gsl_vector_set(u, i, 1./(2.*r + 1.) *

92 (r * (gsl_vector_get(u00, i-1) +

93 gsl_vector_get(u00, i+1)) +

94 gsl_vector_get(u0, i) +

95 + ht * f(t, gsl_vector_get(x, i))));

97 gsl_vector_set(diff, i, gsl_vector_get(u, i) -

98 gsl_vector_get(u00, i));

99 }

100

101 #pragma omp single

102 {

103 norm = gsl_blas_dnrm2(diff);

104 printf("\t%d: norm=%.2e\n", it, norm);

105

106 gsl_vector_memcpy(u00, u);

107 }

108

109 if (norm < tol) {

110 break;

111 }

112

113 }

114

115 #pragma omp single

116 {

117 sprintf(str, "results/u_%06d.dat", k+1);

118 printf("t = %g, [%s]\n", t, str);

119 fu = fopen(str, "wb");

120 gsl_vector_fwrite(fu, u0);

121 fclose(fu);

122

123 // u0 = u

124 gsl_vector_memcpy(u0, u);

125 }

126

127 }

128

129 }

130

131 gsl_vector_free(x);

132 gsl_vector_free(u);

133 gsl_vector_free(u0);

134 gsl_vector_free(u00);

135 gsl_vector_free(diff);

136

137 return 0;

138}

2.6.3 Exercícios

Em construção

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen

| | | |