Colabore! Saiba mais screen_rotation

Ao navegar por este site , você concorda com a Política de Uso de Dados.

| | | | |

Minicurso de C/C++ para Matemática 4 Elementos da Programação Estruturada 6 Elementos da Orientação-a-Objetos

5 Elementos da Computação Matricial

GSL (GNU Scientific Library) é uma biblioteca de métodos numéricos para C/C++. É um software livre sobre a GNU General Public License e está disponível em

https://www.gnu.org/software/gsl/.

A biblioteca fornece uma grande número de rotinas matemáticas para várias áreas da análise numérica. Aqui, vamos nos concentrar em uma rápida introdução à computação matricial.

5.1 Vetores

A alocação de um vetor no GSL segue o estilo de malloc (alocação de memória) e free (liberação de memória). O suporte a vetores requer a importação da biblioteca gsl_vector.h.

Exemplo 5.1.

No seguinte código, alocamos e imprimimos o seguinte vetor

\boldsymbol{v}=(\sqrt{2},1,3.5,\pi).

(18)

Código 9: vector.cc

⬇

1#include <stdio.h>

3// GSL const e funs matemáticas

4#include <gsl/gsl_math.h>

5// GSL vetores

6#include <gsl/gsl_vector.h>

8int main() {

9 // alocação de memória

10 gsl_vector *v = gsl_vector_alloc(4);

12 // atribuição

13 gsl_vector_set(v, 0, sqrt(2.));

14 gsl_vector_set(v, 1, 1.);

15 gsl_vector_set(v, 2, 3.5);

16 gsl_vector_set(v, 3, M_PI);

18 // acesso e impressão

19 for (int i=0; i<4; ++i) {

20 printf("v_%d = %g\n", i, gsl_vector_get(v, i));

21 }

23 // liberação de memória

24 gsl_vector_free(v);

25}

A compilação desse código requer a lincagem com a biblioteca GSL:

⬇

1$ gcc vetor.cc -lgsl -lgslcblas -lm

Observação 5.1.

(Inicialização.) Alternativamente, a alocação com o método

⬇

1gsl_vector *gsl_vector_calloc(size_t n)

cria um vetor e inicializada todos os seus elementos como zero. Outros métodos de inicialização estão disponíveis, consulte GSL Docs: Initializing vector elements.

5.1.1 Operações com Vetores

Operações básicas envolvendo vetores do GSL estão disponíveis com os seguintes métodos¹⁷¹⁷endnote: ¹⁷Mais detalhes, consulte GNU Docs: Vector operations.:

•

int gsl_vector_add(gsl_vector *a, const gsl_vector *b)

Computa a adição vetorial a + b e o resultado é armazenado no vetor a.
•

int gsl_vector_sub(gsl_vector *a, const gsl_vector *b)

Computa a subtração vetorial a - b e o resultado é armazenado no vetor a.
•

int gsl_vector_mul(gsl_vector *a, const gsl_vector *b)

Computa a multiplicação elemento-a-elemento a*b e armazena o resultado no vetor a.
•

int gsl_vector_div(gsl_vector *a, const gsl_vector *b)

Computa a divisão elemento-a-elemento a/b e armazena o resultado no vetor a.
•

int gsl_vector_scale(gsl_vector *a, const double x)

Computa a multiplicação por escalar x*a e armazena o resultado no vetor a.
•

int gsl_vector_add_constant(gsl_vector *a, const double x)

Re-computa o vetor a somando o escalar x a cada um de seus elementos.
•

double gsl_vector_sum(const gsl_vector *a)

Retorna a soma dos elementos do vetor a.

Observação 5.2.

(Suporte BLAS.) O GSL também fornece suporte a biblioteca BLAS (Basic Linear Algebra Subprograms) pelo gsl_blas.h. Operações vetoriais estão disponíveis no nível 1 da biblioteca. Para mais informações e acesso à documentação sobre os métodos disponíveis, consulte

https://www.gnu.org/software/gsl/doc/html/blas.html.

Exemplo 5.2.

No código abaixo, computamos $\boldsymbol{w}=\alpha\boldsymbol{u}+\boldsymbol{v}$ para o escalar $\alpha=2$ e os vetores

		$\displaystyle\boldsymbol{u}=(1,-2,0.5),$		(19)
		$\displaystyle\boldsymbol{v}=(2,1,-1.5).$		(19)

Código 10: axpy.cc

⬇

1#include <stdio.h>

3// GSL vetores

4#include <gsl/gsl_vector.h>

5// GSL BLAS

6#include <gsl/gsl_blas.h>

8int main() {

10 // alpha

11 double alpha = 2.;

13 // u

14 gsl_vector *u = gsl_vector_alloc(3);

15 gsl_vector_set(u, 0, 1.);

16 gsl_vector_set(u, 1, -2.);

17 gsl_vector_set(u, 2, 0.5);

19 // v

20 gsl_vector *v = gsl_vector_alloc(3);

21 gsl_vector_set(v, 0, 2.);

22 gsl_vector_set(v, 1, 1.);

23 gsl_vector_set(v, 2, -1.5);

25 // w = alpha*u + v

26 // alloc w

27 gsl_vector *w = gsl_vector_alloc(3);

28 // copy w = v

29 gsl_vector_memcpy(w, v);

30 // w = alpha*u + w

31 gsl_blas_daxpy(alpha, u, w);

33 // imprime

34 for (int i=0; i<3; ++i) {

35 printf("w_%d = %g\n", i, gsl_vector_get(w, i));

36 }

38 // liberação de memória

39 gsl_vector_free(v);

40}

Exercício 5.1.1.

Faça um código para computar o produto escalar $\boldsymbol{x}\cdot\boldsymbol{y}$ dos vetores

	$\displaystyle\boldsymbol{x}=(1.2,\ln(2),4),$		(20)
	$\displaystyle\boldsymbol{y}=(\pi^{2},\sqrt{3},e).$		(21)

a)

Crie sua própria função double dot(const gsl_vector *x, const gsl_vector *y) que recebe os vetores e retorna o produto escalar deles.
b)

Use o método BLAS gsl_blas_dsdot.

Exercício 5.1.2.

Faça um código para computar a norma $L^{2}$ do vetor

\boldsymbol{x}=(1.2,\log_{10}^{2}(2),0.5).

(22)

a)

Crie sua própria função double norm2(const gsl_vector *x) que recebe o vetor e retorna sua norma.
b)

Use o método BLAS gsl_blas_dnrm2.

5.2 Matrizes

Assim como vetores, a alocação de matrizes no GSL segue o estilo malloc e free. O suporte a matrizes requer a importação da biblioteca gsl_matrix.h.

Exemplo 5.3.

No seguinte código, alocamos e imprimimos a matriz

A=\begin{bmatrix}1.5&-1&\pi&2.3\\ \sqrt[3]{25}&2.1&-3.5&3\\ \log 15&0&2.5^{3}&-1.7\end{bmatrix}

(23)

Código 11: matriz.cc

⬇

1#include <stdio.h>

3// GSL

4#include <gsl/gsl_math.h>

5#include <gsl/gsl_matrix.h>

7int main()

9 // alocação

10 gsl_matrix *A = gsl_matrix_alloc(3, 4);

12 // população

13 gsl_matrix_set(A, 0, 0, 1.5);

14 gsl_matrix_set(A, 0, 1, -1.);

15 gsl_matrix_set(A, 0, 2, M_PI);

16 gsl_matrix_set(A, 0, 3, 2.3);

18 gsl_matrix_set(A, 1, 0, cbrt(25.));

19 gsl_matrix_set(A, 1, 1, 2.1);

20 gsl_matrix_set(A, 1, 2, -3.5);

21 gsl_matrix_set(A, 1, 3, 3.);

23 gsl_matrix_set(A, 2, 0, log10(15.));

24 gsl_matrix_set(A, 2, 1, 0.);

25 gsl_matrix_set(A, 2, 2, pow(2.5, 3.));

26 gsl_matrix_set(A, 2, 3, -1.7);

28 // imprime

29 printf("A = \n");

30 for (int i=0; i<A->size1; ++i) {

31 for (int j=0; j<A->size2; ++j)

32 printf("%g ", gsl_matrix_get(A, i, j));

33 printf("\n");

34 }

36 gsl_matrix_free(A);

37 return 0;

38}

Observação 5.3.

(Inicialização de Matrizes.) Matrizes podem ser inicializadas com todos os seus elementos nulos usando

⬇

1gsl_matrix *gsl_matrix_calloc(size_t n1, size_t n2)

Outros métodos de inicialização também estão disponíveis, consulte

https://www.gnu.org/software/gsl/doc/html/vectors.html#initializing-matrix-elements.

5.3 Operações Matriciais

O GSL conta com os seguintes métodos de operações matriciais¹⁸¹⁸endnote: ¹⁸Consulte a lista completa em GSL Docs: Matrix operations.:

•

int gsl_matrix_add(gsl_matrix *a, const gsl_matrix *b)

Computa a adição matricial a + b e armazena o resultado em a.
•

int gsl_matrix_sub(gsl_matrix *a, const gsl_matrix *b)

Computa a subtração matricial a - b e armazena o resultado em a.
•

int gsl_matrix_mul_elements(gsl_matrix *a, const gsl_matrix *b)

Computa a multiplicação elemento-a-elemento a*b e armazena o resultado em a.
•

int gsl_matrix_div_elements(gsl_matrix *a, const gsl_matrix *b)

Computa a divisão elemento-a-elemento a/b e armazena o resultado em a.
•

int gsl_matrix_scale(gsl_matrix *a, const double x)

Computa a multiplicação por escalar x*a e armazena o resultado em a.
•

int gsl_matrix_add_constant(gsl_matrix *a, const double x)

Re-computa a matriz a somando x a cada um de seus elementos.

Observação 5.4.

(Operações Matrix-Vetor e Matriz-Matriz.) Na GSL, operações matrix-vetor estão disponíveis no suporte BLAS de nível 2. Já, operações matriz-matriz, no suporte BLAS de nível 3. Consulte a lista completa de métodos em

https://www.gnu.org/software/gsl/doc/html/blas.html#blas-support.

Exemplo 5.4.

O seguinte código verifica se $\boldsymbol{x}=(-1,1,-2)$ é solução do sistema linear

		$\displaystyle x_{1}-x_{2}+2x_{3}=-6$		(24)
		$\displaystyle 2x_{1}+x_{2}-x_{3}=1$
		$\displaystyle-x_{1}+x_{2}+x_{3}=0$

Código 12: sisLin.cc

⬇

1#include <stdio.h>

3// GSL

4#include <gsl/gsl_math.h>

5#include <gsl/gsl_vector.h>

6#include <gsl/gsl_matrix.h>

7#include <gsl/gsl_blas.h>

9int main()

10{

11 // matriz dos coefs

12 gsl_matrix *A = gsl_matrix_alloc(3, 3);

14 gsl_matrix_set(A, 0, 0, 1.);

15 gsl_matrix_set(A, 0, 1, -1.);

16 gsl_matrix_set(A, 0, 2, 2.);

18 gsl_matrix_set(A, 1, 0, 2.);

19 gsl_matrix_set(A, 1, 1, 1.);

20 gsl_matrix_set(A, 1, 2, -1.);

22 gsl_matrix_set(A, 2, 0, -1.);

23 gsl_matrix_set(A, 2, 1, 1.);

24 gsl_matrix_set(A, 2, 2, 1.);

26 // vetor dos termos consts

27 gsl_vector *b = gsl_vector_alloc(3);

29 gsl_vector_set(b, 0, -6.);

30 gsl_vector_set(b, 1, 1.);

31 gsl_vector_set(b, 2, 0.);

33 // vetor solução ?

34 gsl_vector *x = gsl_vector_alloc(3);

36 gsl_vector_set(x, 0, -1.);

37 gsl_vector_set(x, 1, 1.);

38 gsl_vector_set(x, 2, -2.);

40 // verificação

41 // y = Ax

42 gsl_vector *y = gsl_vector_alloc(3);

43 gsl_blas_dgemv(CblasNoTrans, 1., A, x, 0., y);

45 // y - b

46 gsl_vector_sub(y, b);

48 if (gsl_blas_dnrm2(y) < 1e-14)

49 printf("x é solução do sistema.\n");

50 else

51 printf("x não é solução do sistema.\n");

53 gsl_matrix_free(A);

54 gsl_vector_free(b);

55 gsl_vector_free(x);

56 gsl_vector_free(y);

58 return 0;

59}

Exercício 5.3.1.

Crie uma função para computar a norma de Frobenius de uma matriz $A$ . Teste seu código com a matriz

A=\begin{bmatrix}-4&-3&-2\\ -1&0&1\\ 2&3&4\end{bmatrix}.

(25)

Exercício 5.3.2.

Faça um código para verificar se a matriz

A=\begin{bmatrix}1&0&0\\ 1&3&1\\ 1&2&0\end{bmatrix}

(26)

é inversa de

B=\begin{bmatrix}1&0&0\\ -0.5&0&0.5\\ 0.5&1&-1.5\end{bmatrix}.

(27)

5.4 Sistemas Lineares

A GSL tem suporte à álgebra linear pelo módulo gsl_linalg.h. Métodos para a decomposição LU, QR, Cholesky entre tantos outros métodos estão disponíveis. Consulte a lista completa em

https://www.gnu.org/software/gsl/doc/html/linalg.html#linear-algebra.

Exemplo 5.5.

No seguinte código, computamos a solução do sistema

		$\displaystyle x_{1}-x_{2}+2x_{3}=-6$		(28)
		$\displaystyle 2x_{1}+x_{2}-x_{3}=1$
		$\displaystyle-x_{1}+x_{2}+x_{3}=0$

pelo método LU.

Código 13: lu.cc

⬇

1#include <stdio.h>

3// GSL

4#include <gsl/gsl_vector.h>

5#include <gsl/gsl_matrix.h>

6#include <gsl/gsl_linalg.h>

8int main()

10 // matriz dos coefs

11 gsl_matrix *A = gsl_matrix_alloc(3, 3);

13 gsl_matrix_set(A, 0, 0, 1.);

14 gsl_matrix_set(A, 0, 1, -1.);

15 gsl_matrix_set(A, 0, 2, 2.);

17 gsl_matrix_set(A, 1, 0, 2.);

18 gsl_matrix_set(A, 1, 1, 1.);

19 gsl_matrix_set(A, 1, 2, -1.);

21 gsl_matrix_set(A, 2, 0, -1.);

22 gsl_matrix_set(A, 2, 1, 1.);

23 gsl_matrix_set(A, 2, 2, 1.);

25 // vetor dos termos consts

26 gsl_vector *b = gsl_vector_alloc(3);

28 gsl_vector_set(b, 0, -6.);

29 gsl_vector_set(b, 1, 1.);

30 gsl_vector_set(b, 2, 0.);

32 // decomposição LU

33 // PA = LU

34 gsl_permutation *p = gsl_permutation_alloc(3);

35 int signum;

36 gsl_linalg_LU_decomp(A, p, &signum);

38 // solução

39 gsl_vector *x = gsl_vector_alloc(3);

40 gsl_linalg_LU_solve(A, p, b, x);

42 // imprime a solução

43 for (int i=0; i<3; ++i)

44 printf("x_%d = %g\n", i, gsl_vector_get(x, i));

46 gsl_matrix_free(A);

47 gsl_vector_free(b);

48 gsl_permutation_free(p);

49 gsl_vector_free(x);

51 return 0;

52}

Exercício 5.4.1.

Crie sua própria função para a computação da solução de sistemas triangulares inferiores. Verifique seu código para o sistema

		$\displaystyle-x_{1}=2$		(29)
		$\displaystyle-3x_{1}+2x_{2}=-8$
		$\displaystyle-x_{1}+x_{2}-x_{3}=0$

Exercício 5.4.2.

Crie sua própria função para a computação da solução de um sistemas triangulares superiores. Verifique seu código para o sistema

		$\displaystyle 2x_{1}-x_{2}+2x_{3}=7$		(30)
		$\displaystyle 2x_{2}-x_{3}=-3$
		$\displaystyle 3x_{3}=3$

Exercício 5.4.3.

Faça um código para resolver o sistema linear

		$\displaystyle x_{1}-x_{2}+2x_{3}=-6$		(31)
		$\displaystyle 2x_{1}+x_{2}-x_{3}=1$
		$\displaystyle-x_{1}+x_{2}+x_{3}=0$

na sua forma matricial $A\boldsymbol{x}=\boldsymbol{b}$ .

a)

Use int gsl_linalg_LU_decomp(gsl_matrix *A, gsl_permutation *p, int *signum) para computar a decomposição $A=LU$ , onde $A$ é a matriz de coeficientes do sistema.
b)

Use sua função criada no Exercício 5.4.1 para resolver $L\boldsymbol{y}=\boldsymbol{b}$ , onde $\boldsymbol{b}$ o vetor dos termos constantes do sistema.
c)

Use sua função criada no Exercício 5.4.2 para resolver $U\boldsymbol{x}=\boldsymbol{y}$ .

Exercício 5.4.4.

Faça um código para computar a solução do sistema

		$\displaystyle 2x_{1}-x_{2}=0$		(32)
		$\displaystyle x_{i-1}-6x_{i}+4x_{i+1}=\operatorname{sen}(\pi i/[2(n-1)])$
		$\displaystyle x_{n-1}+x_{n}=1$

para $i=2,3,\dotsc,n-1$ , $n\geq 3$ . Dica: use o método gsl_linalg_solve_tridiag.

Envie seu comentário

As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a Política de Use de Dados para mais informações. Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

| | | |