Política de uso de dados

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

| | | |

Minicurso de C++ para Matemática

Minicurso de C++ para Matemática 3 Elementos da programação estruturada 5 Elementos da orientação-a-objetos

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

4 Elementos da computação matricial

Vamos explorar a Eigen, biblioteca para computação matricial e tratamento numérico de dados. um software livre sob licença MPL2 e está disponível em https://eigen.tuxfamily.org/. A biblioteca é header-only, ou seja, não requer compilação e instalação. Para incluí-la em um código C++, basta adicionar o módulo eigen3/Eigen/Eigen no cabeçalho.

⬇

1#include <eigen3/Eigen/Eigen>

4.1 Arranjos, matrizes e vetores

Eigen fornece objetos para arranjos de elementos (1D e 2D), matrizes e vetores (matriz linha ou coluna). Os padrões das classes são:

•

Eigen::ArrayXd

Objeto da classe dos arranjos de tamanho dinâmico e de elementos do tipo double. Objeto para organização e processamento de dados com operações elemento-a-elemento.

⬇

1typedef Array<double, Dynamic, 1> ArrayXd;

Código 10: Ex_ArrayXd.cpp

⬇

1// Array 1D de 3 elementos double

2Eigen::ArrayXd x(3);

3// Inicialização

4x << M_PI, M_E, 1.0;

5// atribuição de elemento

6x(2) = 2.0;

7// Impressão

8std::cout << "x = \n" << x << std::endl;
•

Eigen::ArrayXXd

Objeto da classe dos arranjos bidimensionais de tamanho dinâmico e de elementos do tipo double.

⬇

1typedef Array<double, Dynamic, Dynamic> ArrayXXd;

Código 11: Ex_ArrayXXd.cpp

⬇

1// Array 2D 2x3

2Eigen::ArrayXXd a(2, 3);

3// Inicialização

4a << 0.0, 0.1, 0.2,

5 1.0, 1.1, 1.2;

6// atribuição de elemento

7a(1,0) = 2.0;

8// Impressão

9std::cout << "a = \n" << a << std::endl;
•

Eigen::VectorXd

Objeto da classe dos vetores coluna (matriz n x 1) de tamanho dinâmico e de elementos do tipo double. Objeto para modelagem computacional de vetores colunas com operações matriciais.

⬇

1typedef Matrix<double, Dynamic, 1> VectorXd;

Código 12: Ex_VectorXd.cpp

⬇

1// Vetor coluna

2Eigen::VectorXd v {{1.0, 2.0, 3.0}};

3// dimensões

4std::cout << "dimensões de v: "

5 << v.rows() << "x"

6 << v.cols() << std::endl;

7// imprime

8std::cout << "v =\n" << v << std::endl;
•

Eigen::RowVectorXd

Objetos da classe dos vetores linha de tamanho dinâmico e de elementos do tipo double.

⬇

1typedef Matrix<double, 1, Dynamic> RowVectorXd;

Código 13: Ex_RowVectorXd.cpp

⬇

1// Vetor linha

2Eigen::RowVectorXd w(3);

3w << 1.0, 2.0, 3.0;

4// dimensões

5std::cout << "Tamanho de w: "

6 << w.size() << std::endl;

7// imprime

8std::cout << "w =\n" << w << std::endl;
•

Eigen::MatrixXd

Objetos da classe das matrizes de tamanho dinâmico e de elementos do tipo double. Objeto para modelagem computacional com operações matriciais.

⬇

1typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

Código 14: Ex_MatrixXd.cpp

⬇

1// Matriz 3x2

2Eigen::MatrixXd m(3, 2);

3// Inicialização

4m << 1.0, 2.0,

5 3.0, 4.0,

6 5.0, 6.0;

7m(0,0) += m(2,1);

8// Impressão

9std::cout << "m = \n" << m << std::endl;

Observação 4.1.1.(Outros tipos de elementos)

Os elementos das matrizes e arranjos podem ser de outros tipos, como int ou float. Por conveniência, temos typedefs para os tipos mais comuns. Por exemplo, para objetos com elementos int temos ArrayXi, ArrayXXi, VectorXi, RowVectorXi e MatrixXi.

Exercício 4.1.1.

Faça um código que aloque os vetores

	$\displaystyle\vec{u}=(1.2,-3.1,4),$		(13)
	$\displaystyle\vec{v}=(\pi,\sqrt{2},e^{-2})$		(14)

e compute o vetor $\vec{w}=(w_{1},w_{2},w_{3})$ de elementos $w_{i}=u_{i}v_{i}$ , $i=1,2,3$ .

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::VectorXd u(3);

8 u << 1.2, -3.1, 4.0;

10 Eigen::VectorXd v(3);

11 v << M_PI, M_SQRT2, 1.0/(M_E*M_E);

13 Eigen::VectorXd w(u.size());

15 for (int i=0; i<u.size(); i++)

16 {

17 w(i) = u(i) * v(i);

18 }

20 std::cout << "w =\n"

21 << w << std::endl;

23 return 0;

24}

Exercício 4.1.2.

Crie um código que aloque a matriz e o vetor

	$\displaystyle A=\begin{bmatrix}-1&\frac{1}{3}\\[6.0pt]2&\sqrt{2}/2\\[6.0pt]e^{-1}&-3\end{bmatrix},$		(15)
	$\displaystyle\vec{v}=[-5,\pi/2,\operatorname{sen}(\pi/3)].$		(16)

e compute o vetor $\vec{w}=[w_{1},w_{2}]$ de elementos

w_{j}=\sum_{i=1}^{3}v_{i}A_{ij}.

(17)

para $j=1,2$ .

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::MatrixXd A {{-1.0, 1.0/3,},

8 { 2.0, M_SQRT1_2},

9 { 1.0/M_E, -3.0}};

11 Eigen::RowVectorXd v(3);

12 v << -5.0, M_PI_2, sin(M_PI/3);

14 Eigen::RowVectorXd w(A.cols());

16 for (int j=0; j<A.cols(); j++) {

17 w(j) = 0.0;

18 for (int i=0; i<A.rows(); i++) {

19 w(j) += v(j) * A(i,j);

20 }

21 }

23 std::cout << "w =\n"

24 << w << std::endl;

26 return 0;

27}

Exercício 4.1.3.

Implemente uma função que computa o determinante de matrizes reais A $3\times 3$ . Compare o resultado com o método A.determinant() da Eigen.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

5double det(Eigen::Matrix3d &A)

7 // Determinante de uma matriz 3x3

8 return A(0, 0) * (A(1, 1) * A(2, 2) - \

9 A(1, 2) * A(2, 1)) -\

10 A(0, 1) * (A(1, 0) * A(2, 2) - \

11 A(1, 2) * A(2, 0)) +\

12 A(0, 2) * (A(1, 0) * A(2, 1) - \

13 A(1, 1) * A(2, 0));

14}

16int main()

17{

18 Eigen::Matrix3d A;

19 A << -1.0, -2.0, -1.0,

20 2.0, 1.0, 0.0,

21 1.0, -1.0, -2.0;

23 std::cout << "A =\n" << std::endl

24 << A << std::endl;

25 std::cout << "det(A) = "

26 << det(A) << std::endl;

28 return 0;

29}

4.2 Blocos e fatiamentos

4.2.1 Blocos

Eigen permite acessar blocos de matrizes e vetores com o método block(i, j, n, m), onde i e j são as coordenadas iniciais do bloco e n e m são as dimensões do bloco. O método block retorna um objeto do tipo Eigen::Block que pode ser manipulado como uma matriz ou vetor.

⬇

1Eigen::MatrixXd A(4,4);

2A << 1, 2, 3, 4,

3 5, 6, 7, 8,

4 9, 10, 11, 12,

5 13, 14, 15, 16;

7Eigen::Block block = A.block(1,0,3,2);

8std::cout << "A[1:3, 0:2] =\n"

9 << block << std::endl;

11block(1,1) = 99;

12std::cout << "A =\n"

13 << A << std::endl;

A[1:3, 0:2] =

5 6

9 10

13 14

A =

1 2 3 4

5 6 7 8

9 99 11 12

13 14 15 16

4.2.2 Linhas e colunas

Eigen permite acessar linhas e colunas de matrizes com o método row(i) e col(j), respectivamente. O método row e col retorna um objeto do tipo Eigen::Matrix que pode ser manipulado como um vetor ou matriz.

⬇

1Eigen::MatrixXd A(4,4);

2A << 1, 2, 3, 4,

3 5, 6, 7, 8,

4 9, 10, 11, 12,

5 13, 14, 15, 16;

7std::cout << "A[1,:] = "

8 << A.row(1) << std::endl;

10A.col(2) << 0, 0, 0, 0;

11std::cout << "A = " << A << std::endl;

A[1,:] = 5 6 7 8

A = 1 2 0 4

5 6 0 8

9 10 0 12

13 14 0 16

4.2.3 Fatiamentos

O fatiamento de matrizes é convenientemente feito com o ajuda do método Eigen::seq(firstIdx,lastIdx,incr), uma sequência de inteiros de firstIdx até lastIdx com incremento incr. Por exemplo, Eigen::seq(0,10,2) retorna a sequência 0, 2, 4, 6, 8, 10. Estudemos o seguinte código.

Código 15: Ex_slice.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd m(3,4);

7 m << 1, 2, 3, -1,

8 4, 5, 6, -2,

9 7, 8, 9, -3;

10 std::cout << "m =" << std::endl

11 << m << std::endl;

13 std::cout << "m[1:2,0:3:2] =" << std::endl

14 << m(Eigen::seq(1,2), Eigen::seq(0,3,2)) << std::endl;

16 std::cout << "m[:,1:3:2] =" << std::endl

17 << m(Eigen::placeholders::all, Eigen::seq(1,3,2)) << std::endl;

19 Eigen::ArrayXi ind {{0, 1, 3}};

20 std::cout << "m[:,(0,1,3)] =" << std::endl

21 << m(Eigen::placeholders::all, ind) << std::endl;

23 return 0;

24}

m =

1 2 3 -1

4 5 6 -2

7 8 9 -3

m[1:2,0:3:2] =

4 6

7 9

m[:,1:3:2] =

2 -1

5 -2

8 -3

m[:,(0,1,3)] =

1 2 -1

4 5 -2

7 8 -3

Observação 4.2.1.(Mais sobre fatiamentos)

Consulte mais sobre fatiamentos em Eigen::Slicing and Indexing.

Exercício 4.2.1.

Faça um código que aloque a matriz

A=\begin{bmatrix}1&2&3&4\\5&6&7&8\\9&10&11&12\\13&14&15&16\end{bmatrix}

(18)

e imprima os seguintes fatiamentos:

a)

$A[1:3,0:2]$
b)

$A[0:2,1:3:2]$
c)

$A[0:2:2,0:3:3]$
d)

$A[0:3:2,:]$
e)

$A[(0,2,3),1:3]$

Resposta 0.

a) a[1:3,0:2] =

5 6 7

9 10 11

13 14 15

b) a[0:2,1:3:2] =

2 4

6 8

10 12

c) a[0:2:2,0:3:3] =

1 4

9 12

d) a[0:3:2,:] =

1 2 3 4

9 10 11 12

a[(0,2,3),1:3] =

2 3 4

10 11 12

14 15 16

4.3 Operações matriciais

A Eigen oferece operações matriciais como soma, subtração e multiplicação de matrizes pela sobrecarga dos operadores +, - e *, respectivamente.

4.3.1 Adição e subtração

A adição e subtração de matrizes ou vetores é feita com os operadores + e -, respectivamente. A adição e subtração de matrizes e vetores é feita elemento-a-elemento. Para a adição e subtração de matrizes, as dimensões devem ser compatíveis. Para a adição e subtração de vetores, os tamanhos devem ser compatíveis. Se não forem, o programa falha e interrompe sua execução. Estudemos o seguinte código.

Código 16: ExSomaDiferenca.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 // adição de matrizes

7 Eigen::MatrixXd a {{1,2},{3,4}};

8 Eigen::MatrixXd b {{5,6},{7,8}};

10 std::cout << "a + b = \n" <<

11 a + b << std::endl;

13 // subtração de vetores

14 Eigen::VectorXd v {{1,2,3}};

15 Eigen::VectorXd w {{4,5,6}};

17 std::cout << "v + w = \n" <<

18 v - w << std::endl;

21 return 0;

22}

a + b =

6 8

10 12

v + w =

-3

4.3.2 Multiplicação e divisão por escalar

A multiplicação e divisão de matrizes e vetores por escalares é feita com os operadores * e /, respectivamente. Estudemos o seguinte código.

Código 17: ExEscalar.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 double alpha = 0.5;

7 Eigen::MatrixXd a {{1,2},{3,4}};

8 Eigen::VectorXd v {{1,2,3}};

10 std::cout << "2a = \n" <<

11 2*a << std::endl;

13 std::cout << "v/alpha = \n" <<

14 v/alpha << std::endl;

16 return 0;

17}

2a =

2 4

6 8

v/alpha =

4.3.3 Multiplicação de matrizes

A multiplicação de matrizes é feita com o operador *. A multiplicação de matrizes é feita com a regra da álgebra linear. Para a multiplicação de matrizes, as dimensões devem ser compatíveis. Se não forem, o programa falha e interrompe sua execução. Estudemos o seguinte código.

Código 18: ExMultMatrizes.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{1,2,3},

7 {3,4,5}};

8 Eigen::MatrixXd b {{-2,1},

9 {1,-1},

10 {1,1}};

12 std::cout << "a*b = \n" <<

13 a*b << std::endl;

15 Eigen::RowVectorXd u {{1,2,3}};

17 std::cout << "u*b = \n" <<

18 u*b << std::endl;

20 return 0;

21}

a*b =

3 2

3 4

u*b =

3 2

Exercício 4.3.1.

Considere o seguinte sistema de equações lineares

	$\displaystyle 2x_{1}+3x_{2}$	$\displaystyle=-1,$		(19)
	$\displaystyle 3x_{1}-4x_{2}$	$\displaystyle=7.$		(20)

Escreva um código que aloque a matriz dos coeficientes $A$ e o vetor dos termos constantes $\vec{b}$ do sistema. Então, verifique se o vetor $x=(1,-1)$ é solução do sistema.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{2,3},

7 {3,-4}};

8 Eigen::VectorXd b {{-1,7}};

9 Eigen::VectorXd x {{1,-1}};

11 std::cout << "a*x - b = \n"

12 << a*x - b << std::endl;

14 return 0;

15}

Exercício 4.3.2.

Considere o seguinte sistema de equações lineares

$\displaystyle 2x_{1}+3x_{2}+4x_{3}$	$\displaystyle=-1,$	(21)
$\displaystyle 3x_{1}-4x_{2}+5x_{3}$	$\displaystyle=7,$	(22)
$\displaystyle x_{1}+x_{2}-x_{3}$	$\displaystyle=0.$	(23)

Escreva um código que aloque a matriz dos coeficientes $A$ e o vetor dos termos constantes $\vec{b}$ do sistema. Então, aloque a matriz estendida $E=[A|b]$ e compute a equivalente matriz estendida reduzida $E\approx[I|x]$ . Por fim, aloque e imprima o vetor solução $x$ .

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4#define all Eigen::placeholders::all

6int main()

8 // matrix dos coefs

9 Eigen::MatrixXd a {{2,3,4},

10 {3,-4,5},

11 {1,1,-1}};

12 // vetor dos termos constantes

13 Eigen::VectorXd b {{-1,7,0}};

15 // matrix estendida

16 Eigen::MatrixXd ab(a.rows(), a.cols() + 1);

17 ab << a, b;

18 std::cout << "Matriz estendida:\n" << ab << std::endl;

20 // escalonamento

21 ab(1,all) -= (ab(1,0) / ab(0,0)) * ab(0,all);

22 ab(2,all) -= (ab(2,0) / ab(0,0)) * ab(0,all);

23 ab(2,all) -= (ab(2,1) / ab(1,1)) * ab(1,all);

25 std::cout << "Matriz escalonada:\n" << ab << std::endl;

27 // redução

28 ab(1,all) -= (ab(1,2) / ab(2,2)) * ab(2,all);

29 ab(0,all) -= (ab(0,2) / ab(2,2)) * ab(2,all);

30 ab(0,all) -= (ab(0,1) / ab(1,1)) * ab(1,all);

31 ab(1,all) /= ab(1,1);

32 ab(2,all) /= ab(2,2);

33 ab(0,all) /= ab(0,0);

35 std::cout << "Matriz reduzida:\n" << ab << std::endl;

37 // solução

38 Eigen::VectorXd x(ab.rows());

39 x << ab(all, ab.cols()-1);

41 std::cout << "Solução x = \n" << x << std::endl;

42 return 0;

43}

4.4 Operações elemento-a-elemento

Operações aritméticas elemento-a-elemento podem ser feitas com objetos da classe Eigen::Array (Eigen::ArrayXd e Eigen::ArrayXXd, por exemplo). Para matrizes é necessário usar o método array() para converter a matriz em um objeto do tipo Eigen::Array. Estudemos o seguinte código.

Código 19: ExOpElem.cpp

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 // mult. elemento-a-elemento

8 Eigen::ArrayXd arr {{1,2,3}};

9 Eigen::VectorXd vec {{4,5,6}};

11 std::cout << "arr * vec.array() =\n"

12 << arr * vec.array() << std::endl;

14 // divisão elemento-a-elemento

15 Eigen::MatrixXd a {{3, 6}, {12, 21}};

16 Eigen::MatrixXd b {{3, 2}, {2, 7}};

18 std::cout << "a.array() * b.array() =\n"

19 << a.array() / b.array() << std::endl;

21 // funções elementares

22 arr.resize(5);

23 arr << 0.0 , M_PI/3, M_PI_4,

24 M_PI/3, M_PI_2;

25 std::cout << "sin(arr) =\n"

26 << arr.sin() << std::endl;

28 std::cout << "b**2 =\n"

29 << b.array().pow(2) << std::endl;

31 return 0;

32}

arr * vec.array() =

a.array() * b.array() =

1 3

6 3

sin(arr) =

0.866025

0.707107

0.866025

b**2 =

9 4

4 49

Observação 4.4.1.(Mais sobre operações elemento-a-elemento)

Consulte mais sobre operações elemento-a-elemento em Eigen::The Array class and coefficient-wise operations.

Exercício 4.4.1.

Faça um código para computar os valores da função cosseno aplicada ao vetor

\boldsymbol{\theta}=(0,30^{\circ},45^{\circ},60^{\circ},90^{\circ}).

(24)

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::VectorXd theta {{0.0,

8 M_PI/6,

9 M_PI/4,

10 M_PI/3,

11 M_PI/2}};

12 std::cout << "cos(theta) = " << std::endl

13 << theta.array().cos() << std::endl;

15 return 0;

16}

Exercício 4.4.2.

Considere que os seguintes objetos estão declarados declarados:

⬇

1Eigen::ArrayXd arr {{1,2,3}};

2Eigen::VectorXd vec {{4,5,6}};

3Eigen::MatrixXd a {{3, 6, 4},

4 {2, 5, 3}};

Antes de implementar para verificar, forneça o resultado, quando definido, de cada uma das seguintes operações:

a)

arr - vec.array();
b)

a.array() * vec;
c)

a * vec;
d)

a.row(0).array() * arr;
e)

a.col(1) * vec;
f)

a.col(1).array() * arr;

Resposta 0.

a) $(-3,-3,-3)$ ; b) não definido; c) $(66,51)$ ; d) não definido.; e) $51$ ; f) não definido.

4.5 Inicialização

A biblioteca Eigen conta com métodos para inicialização de matrizes e arranjos.

•

Zero()

⬇

1arr = Eigen::ArrayXd::Zero(3)

⬇

arr =

0

0

0

⬇

1v = Eigen::RowVectorXd::Zero(2)

⬇

v = 0 0

⬇

1a = Eigen::MatrixXd::Zero(2,3)

⬇

a =

0 0 0

0 0 0
•

Ones()

⬇

1arr = Eigen::ArrayXd::Ones(2)

⬇

arr =

1

1
•

Constant()

⬇

1a = Eigen::MatrixXd::Constant(2,3,0.5)

⬇

a =

0.5 0.5 0.5

0.5 0.5 0.5
•

Random()

⬇

1a = Eigen::MatrixXd::Random(3,2)

⬇

a =

0.680375 0.59688

-0.211234 0.823295

0.566198 -0.604897
•

Identity()

⬇

1a = Eigen::MatrixXd::Identity(3,3)

⬇

a =

1 0 0

0 1 0

0 0 1
•

LinSpaced()

⬇

1arr = Eigen::RowVectorXd::LinSpaced(5, 0, 1)

⬇

arr =

0 0.25 0.5 0.75 1

Exercício 4.5.1.

Considere as seguintes partições uniformes de intervalos

	$\displaystyle x_{[0,1]}=\{0=x_{0}<x_{1}<\cdots<x_{n_{x}+1}=1\},$		(25)
	$\displaystyle y_{[1,2]}=\{1=y_{0}<y_{1}<\cdots<y_{n_{y}+1}=2\},$		(26)

onde $n_{x}=5$ e $n_{y}=3$ . Implemente uma função X, Y = meshgrid(x, y), em que X e Y são os arranjos $n_{x}\times n_{y}$ que contém as coordenadas dos pontos do produto cartesiano $x_{[0,1]}\times y_{[1,2]}$ . Isto é, a grade de pontos $(X_{i,j},Y_{i,j})=(x_{i},y_{j})$ , $i=0,1,\ldots,n_{x}$ e $j=0,1,\ldots,n_{y}$ .

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4void meshgrid(const Eigen::ArrayXd& x, const Eigen::ArrayXd& y,

5 Eigen::ArrayXXd& X, Eigen::ArrayXXd& Y)

7 X = x.replicate(1, y.size());

8 Y = y.transpose().replicate(x.size(), 1);

11int main()

12{

13 Eigen::ArrayXd x = Eigen::ArrayXd::LinSpaced(5, 0, 1);

14 Eigen::ArrayXd y = Eigen::ArrayXd::LinSpaced(3, 1, 2);

16 Eigen::ArrayXXd X;

17 Eigen::ArrayXXd Y;

18 meshgrid(x, y, X, Y);

20 std::cout << "X:\n" << X << "\n";

21 std::cout << "Y:\n" << Y << "\n";

23 return 0;

24}

Exercício 4.5.2.

Implemente uma função para computar a inverse de uma matriz não-singular pelo método de eliminação gaussiana. Verifique seu código com o método A.inverse() da Eigen.

Resposta 0.

⬇

1Eigen::MatrixXd inverse(const Eigen::MatrixXd& A)

3 // extendida

4 Eigen::MatrixXd E(A.rows(), 2 * A.cols());

5 E << A, Eigen::MatrixXd::Identity(A.rows(), A.cols());

7 // escalona

8 for (int j = 0; j < A.cols()-1; ++j) {

9 for (int i = j+1; i < A.rows(); ++i) {

10 E.row(i) -= E(i,j)/E(j,j) * E.row(j);

11 }

12 }

13 // reduzida

14 for (int j = A.cols()-1; j > 0; --j) {

15 for (int i = j-1; i >= 0; --i) {

16 E.row(i) -= E(i,j)/E(j,j) * E.row(j);

17 }

18 }

19 // normaliza

20 for (int j = 0; j < A.cols(); ++j) {

21 E.row(j) /= E(j,j);

22 }

24 return E.rightCols(A.cols());

25}

4.6 Reduções

Reduções são operações que transformam uma matriz ou um arranjo em um escalar.

•

sum()

⬇

1Eigen::MatrixXd a {{1,2,3},

2 {4,5,6}};

3double s = a.sum();

4std::cout << "s = " << s << std::endl;

⬇

s = 21
•

prod()

⬇

1double p = a.prod();

2std::cout << "p = " << p << std::endl;

⬇

p = 720
•

mean()

⬇

1double m = a.mean();

2std::cout << "m = " << m << std::endl;

⬇

m = 3.5
•

minCoeff()/maxCoeff()

⬇

1double min = a.minCoeff();

2std::cout << "min = " << min << std::endl;

3double max = a.maxCoeff();

4std::cout << "max = " << max << std::endl;

⬇

min = 1

max = 6
•

trace()

⬇

1Eigen::MatrixXd b {{1,3},

2 {4,6}};

3

4double tr = b.trace();

5std::cout << "tr = " << tr << std::endl;

⬇

tr = 7
•

norm()

⬇

1Eigen::VectorXd v {{1,2,3}};

2Eigen::MatrixXd m {{1,2},

3 {3,4}};

4// norma l2

5double v_norm = v.norm();

6std::cout << "v_norm = " << v_norm << std::endl;

7// norma L2 (Frobenius)

8double m_norm = m.norm();

9std::cout << "m_norm = " << m_norm << std::endl;

⬇

n = 3.74166

m_norm = 5.47723

Exercício 4.6.1.

Faça um código que aloque os seguintes vetores

	$\displaystyle\vec{u}=(1,2,3),$		(28)
	$\displaystyle\vec{v}=(4,5,6),$		(29)

e compute a o produto interno (ou escalar) $u\cdot v$ . Verifique seu resultado com o método u.dot(v) da Eigen.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::VectorXd u {{1,2,3}};

7 Eigen::VectorXd v {{4,5,6}};

9 Eigen::VectorXd w(u.size());

10 w << u.array() * v.array();

11 std::cout << "u.v = " << w.sum() << std::endl;

12 std::cout << "u.dot(v) = " << u.dot(v) << std::endl;

14 return 0;

15}

Exercício 4.6.2.

Implemente uma função para computar o erro médio quadrático entre dois vetores u e v de mesmo tamanho. O erro médio quadrático é dado por

\text{EMQ}=\frac{1}{n}\sum_{i=1}^{n}(u_{i}-v_{i})^{2}.

(30)

Resposta 0.

⬇

1double EMQ(const Eigen::VectorXd& u, const Eigen::VectorXd& v)

3 return ((u - v).array().pow(2)).sum() / u.size();

4.7 Sistemas lineares

A biblioteca Eigen contém a implementação de métodos de decomposição (LU, QR, SVD) que podem ser usadas para computar a solução de sistemas lineares $Ax=b$ . Por exemplo, no seguinte código computa a solução do sistema

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

usando a o método da decomposição LU com pivotamento parcial.

Código 20: ExSistLin.cpp

⬇

1Eigen::MatrixXd m {{2,-1},

2 {1,3}};

3Eigen::VectorXd b {{3,-2}};

5std::cout << "x =\n"

6 << m.lu().solve(b) << std::endl;

x =

-1

Observação 4.7.1.(Mais sobre sistemas lineares)

Consulte mais sobre sistemas lineares em Eigen::Linear Algebra.

Exercício 4.7.1.

Considere o seguinte sistema de equações lineares

$\displaystyle 2x_{1}+3x_{2}+4x_{3}$	$\displaystyle=-1,$	(33)
$\displaystyle 3x_{1}-4x_{2}+5x_{3}$	$\displaystyle=7,$	(34)
$\displaystyle x_{1}+x_{2}-x_{3}$	$\displaystyle=0.$	(35)

Escreva um código que aloque a matriz dos coeficientes $A$ e o vetor dos termos constantes $\vec{b}$ do sistema. Então, compute a solução do sistema usando o método LU com pivotamento parcial.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{2,3,4},

7 {3,-4,5},

8 {1,1,-1}};

9 Eigen::VectorXd b {{-1,7,0}};

11 Eigen::VectorXd x = a.lu().solve(b);

12 std::cout << "x =\n" << x << std::endl;

14 return 0;

15}

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.

Minicurso de C++ para Matemática

Minicurso de C++ para Matemática 3 Elementos da programação estruturada 5 Elementos da orientação-a-objetos

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

4 Elementos da computação matricial

⬇

1#include <eigen3/Eigen/Eigen>

4.1 Arranjos, matrizes e vetores

Eigen fornece objetos para arranjos de elementos (1D e 2D), matrizes e vetores (matriz linha ou coluna). Os padrões das classes são:

•

Eigen::ArrayXd

Objeto da classe dos arranjos de tamanho dinâmico e de elementos do tipo double. Objeto para organização e processamento de dados com operações elemento-a-elemento.

⬇

1typedef Array<double, Dynamic, 1> ArrayXd;

Código 10: Ex_ArrayXd.cpp

⬇

1// Array 1D de 3 elementos double

2Eigen::ArrayXd x(3);

3// Inicialização

4x << M_PI, M_E, 1.0;

5// atribuição de elemento

6x(2) = 2.0;

7// Impressão

8std::cout << "x = \n" << x << std::endl;
•

Eigen::ArrayXXd

Objeto da classe dos arranjos bidimensionais de tamanho dinâmico e de elementos do tipo double.

⬇

1typedef Array<double, Dynamic, Dynamic> ArrayXXd;

Código 11: Ex_ArrayXXd.cpp

⬇

1// Array 2D 2x3

2Eigen::ArrayXXd a(2, 3);

3// Inicialização

4a << 0.0, 0.1, 0.2,

5 1.0, 1.1, 1.2;

6// atribuição de elemento

7a(1,0) = 2.0;

8// Impressão

9std::cout << "a = \n" << a << std::endl;
•

Eigen::VectorXd

Objeto da classe dos vetores coluna (matriz n x 1) de tamanho dinâmico e de elementos do tipo double. Objeto para modelagem computacional de vetores colunas com operações matriciais.

⬇

1typedef Matrix<double, Dynamic, 1> VectorXd;

Código 12: Ex_VectorXd.cpp

⬇

1// Vetor coluna

2Eigen::VectorXd v {{1.0, 2.0, 3.0}};

3// dimensões

4std::cout << "dimensões de v: "

5 << v.rows() << "x"

6 << v.cols() << std::endl;

7// imprime

8std::cout << "v =\n" << v << std::endl;
•

Eigen::RowVectorXd

Objetos da classe dos vetores linha de tamanho dinâmico e de elementos do tipo double.

⬇

1typedef Matrix<double, 1, Dynamic> RowVectorXd;

Código 13: Ex_RowVectorXd.cpp

⬇

1// Vetor linha

2Eigen::RowVectorXd w(3);

3w << 1.0, 2.0, 3.0;

4// dimensões

5std::cout << "Tamanho de w: "

6 << w.size() << std::endl;

7// imprime

8std::cout << "w =\n" << w << std::endl;
•

Eigen::MatrixXd

Objetos da classe das matrizes de tamanho dinâmico e de elementos do tipo double. Objeto para modelagem computacional com operações matriciais.

⬇

1typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

Código 14: Ex_MatrixXd.cpp

⬇

1// Matriz 3x2

2Eigen::MatrixXd m(3, 2);

3// Inicialização

4m << 1.0, 2.0,

5 3.0, 4.0,

6 5.0, 6.0;

7m(0,0) += m(2,1);

8// Impressão

9std::cout << "m = \n" << m << std::endl;

Observação 4.1.1.(Outros tipos de elementos)

Exercício 4.1.1.

Faça um código que aloque os vetores

	$\displaystyle\vec{u}=(1.2,-3.1,4),$		(13)
	$\displaystyle\vec{v}=(\pi,\sqrt{2},e^{-2})$		(14)

e compute o vetor $\vec{w}=(w_{1},w_{2},w_{3})$ de elementos $w_{i}=u_{i}v_{i}$ , $i=1,2,3$ .

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::VectorXd u(3);

8 u << 1.2, -3.1, 4.0;

10 Eigen::VectorXd v(3);

11 v << M_PI, M_SQRT2, 1.0/(M_E*M_E);

13 Eigen::VectorXd w(u.size());

15 for (int i=0; i<u.size(); i++)

16 {

17 w(i) = u(i) * v(i);

18 }

20 std::cout << "w =\n"

21 << w << std::endl;

23 return 0;

24}

Exercício 4.1.2.

Crie um código que aloque a matriz e o vetor

	$\displaystyle A=\begin{bmatrix}-1&\frac{1}{3}\\[6.0pt]2&\sqrt{2}/2\\[6.0pt]e^{-1}&-3\end{bmatrix},$		(15)
	$\displaystyle\vec{v}=[-5,\pi/2,\operatorname{sen}(\pi/3)].$		(16)

e compute o vetor $\vec{w}=[w_{1},w_{2}]$ de elementos

w_{j}=\sum_{i=1}^{3}v_{i}A_{ij}.

(17)

para $j=1,2$ .

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::MatrixXd A {{-1.0, 1.0/3,},

8 { 2.0, M_SQRT1_2},

9 { 1.0/M_E, -3.0}};

11 Eigen::RowVectorXd v(3);

12 v << -5.0, M_PI_2, sin(M_PI/3);

14 Eigen::RowVectorXd w(A.cols());

16 for (int j=0; j<A.cols(); j++) {

17 w(j) = 0.0;

18 for (int i=0; i<A.rows(); i++) {

19 w(j) += v(j) * A(i,j);

20 }

21 }

23 std::cout << "w =\n"

24 << w << std::endl;

26 return 0;

27}

Exercício 4.1.3.

Implemente uma função que computa o determinante de matrizes reais A $3\times 3$ . Compare o resultado com o método A.determinant() da Eigen.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

5double det(Eigen::Matrix3d &A)

7 // Determinante de uma matriz 3x3

8 return A(0, 0) * (A(1, 1) * A(2, 2) - \

9 A(1, 2) * A(2, 1)) -\

10 A(0, 1) * (A(1, 0) * A(2, 2) - \

11 A(1, 2) * A(2, 0)) +\

12 A(0, 2) * (A(1, 0) * A(2, 1) - \

13 A(1, 1) * A(2, 0));

14}

16int main()

17{

18 Eigen::Matrix3d A;

19 A << -1.0, -2.0, -1.0,

20 2.0, 1.0, 0.0,

21 1.0, -1.0, -2.0;

23 std::cout << "A =\n" << std::endl

24 << A << std::endl;

25 std::cout << "det(A) = "

26 << det(A) << std::endl;

28 return 0;

29}

4.2 Blocos e fatiamentos

4.2.1 Blocos

⬇

1Eigen::MatrixXd A(4,4);

2A << 1, 2, 3, 4,

3 5, 6, 7, 8,

4 9, 10, 11, 12,

5 13, 14, 15, 16;

7Eigen::Block block = A.block(1,0,3,2);

8std::cout << "A[1:3, 0:2] =\n"

9 << block << std::endl;

11block(1,1) = 99;

12std::cout << "A =\n"

13 << A << std::endl;

A[1:3, 0:2] =

5 6

9 10

13 14

A =

1 2 3 4

5 6 7 8

9 99 11 12

13 14 15 16

4.2.2 Linhas e colunas

⬇

1Eigen::MatrixXd A(4,4);

2A << 1, 2, 3, 4,

3 5, 6, 7, 8,

4 9, 10, 11, 12,

5 13, 14, 15, 16;

7std::cout << "A[1,:] = "

8 << A.row(1) << std::endl;

10A.col(2) << 0, 0, 0, 0;

11std::cout << "A = " << A << std::endl;

A[1,:] = 5 6 7 8

A = 1 2 0 4

5 6 0 8

9 10 0 12

13 14 0 16

4.2.3 Fatiamentos

Código 15: Ex_slice.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd m(3,4);

7 m << 1, 2, 3, -1,

8 4, 5, 6, -2,

9 7, 8, 9, -3;

10 std::cout << "m =" << std::endl

11 << m << std::endl;

13 std::cout << "m[1:2,0:3:2] =" << std::endl

14 << m(Eigen::seq(1,2), Eigen::seq(0,3,2)) << std::endl;

16 std::cout << "m[:,1:3:2] =" << std::endl

17 << m(Eigen::placeholders::all, Eigen::seq(1,3,2)) << std::endl;

19 Eigen::ArrayXi ind {{0, 1, 3}};

20 std::cout << "m[:,(0,1,3)] =" << std::endl

21 << m(Eigen::placeholders::all, ind) << std::endl;

23 return 0;

24}

m =

1 2 3 -1

4 5 6 -2

7 8 9 -3

m[1:2,0:3:2] =

4 6

7 9

m[:,1:3:2] =

2 -1

5 -2

8 -3

m[:,(0,1,3)] =

1 2 -1

4 5 -2

7 8 -3

Observação 4.2.1.(Mais sobre fatiamentos)

Consulte mais sobre fatiamentos em Eigen::Slicing and Indexing.

Exercício 4.2.1.

Faça um código que aloque a matriz

A=\begin{bmatrix}1&2&3&4\\5&6&7&8\\9&10&11&12\\13&14&15&16\end{bmatrix}

(18)

e imprima os seguintes fatiamentos:

a)

$A[1:3,0:2]$
b)

$A[0:2,1:3:2]$
c)

$A[0:2:2,0:3:3]$
d)

$A[0:3:2,:]$
e)

$A[(0,2,3),1:3]$

Resposta 0.

a) a[1:3,0:2] =

5 6 7

9 10 11

13 14 15

b) a[0:2,1:3:2] =

2 4

6 8

10 12

c) a[0:2:2,0:3:3] =

1 4

9 12

d) a[0:3:2,:] =

1 2 3 4

9 10 11 12

a[(0,2,3),1:3] =

2 3 4

10 11 12

14 15 16

4.3 Operações matriciais

A Eigen oferece operações matriciais como soma, subtração e multiplicação de matrizes pela sobrecarga dos operadores +, - e *, respectivamente.

4.3.1 Adição e subtração

Código 16: ExSomaDiferenca.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 // adição de matrizes

7 Eigen::MatrixXd a {{1,2},{3,4}};

8 Eigen::MatrixXd b {{5,6},{7,8}};

10 std::cout << "a + b = \n" <<

11 a + b << std::endl;

13 // subtração de vetores

14 Eigen::VectorXd v {{1,2,3}};

15 Eigen::VectorXd w {{4,5,6}};

17 std::cout << "v + w = \n" <<

18 v - w << std::endl;

21 return 0;

22}

a + b =

6 8

10 12

v + w =

-3

4.3.2 Multiplicação e divisão por escalar

A multiplicação e divisão de matrizes e vetores por escalares é feita com os operadores * e /, respectivamente. Estudemos o seguinte código.

Código 17: ExEscalar.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 double alpha = 0.5;

7 Eigen::MatrixXd a {{1,2},{3,4}};

8 Eigen::VectorXd v {{1,2,3}};

10 std::cout << "2a = \n" <<

11 2*a << std::endl;

13 std::cout << "v/alpha = \n" <<

14 v/alpha << std::endl;

16 return 0;

17}

2a =

2 4

6 8

v/alpha =

4.3.3 Multiplicação de matrizes

Código 18: ExMultMatrizes.cpp

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{1,2,3},

7 {3,4,5}};

8 Eigen::MatrixXd b {{-2,1},

9 {1,-1},

10 {1,1}};

12 std::cout << "a*b = \n" <<

13 a*b << std::endl;

15 Eigen::RowVectorXd u {{1,2,3}};

17 std::cout << "u*b = \n" <<

18 u*b << std::endl;

20 return 0;

21}

a*b =

3 2

3 4

u*b =

3 2

Exercício 4.3.1.

Considere o seguinte sistema de equações lineares

	$\displaystyle 2x_{1}+3x_{2}$	$\displaystyle=-1,$		(19)
	$\displaystyle 3x_{1}-4x_{2}$	$\displaystyle=7.$		(20)

Escreva um código que aloque a matriz dos coeficientes $A$ e o vetor dos termos constantes $\vec{b}$ do sistema. Então, verifique se o vetor $x=(1,-1)$ é solução do sistema.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{2,3},

7 {3,-4}};

8 Eigen::VectorXd b {{-1,7}};

9 Eigen::VectorXd x {{1,-1}};

11 std::cout << "a*x - b = \n"

12 << a*x - b << std::endl;

14 return 0;

15}

Exercício 4.3.2.

Considere o seguinte sistema de equações lineares

$\displaystyle 2x_{1}+3x_{2}+4x_{3}$	$\displaystyle=-1,$	(21)
$\displaystyle 3x_{1}-4x_{2}+5x_{3}$	$\displaystyle=7,$	(22)
$\displaystyle x_{1}+x_{2}-x_{3}$	$\displaystyle=0.$	(23)

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4#define all Eigen::placeholders::all

6int main()

8 // matrix dos coefs

9 Eigen::MatrixXd a {{2,3,4},

10 {3,-4,5},

11 {1,1,-1}};

12 // vetor dos termos constantes

13 Eigen::VectorXd b {{-1,7,0}};

15 // matrix estendida

16 Eigen::MatrixXd ab(a.rows(), a.cols() + 1);

17 ab << a, b;

18 std::cout << "Matriz estendida:\n" << ab << std::endl;

20 // escalonamento

21 ab(1,all) -= (ab(1,0) / ab(0,0)) * ab(0,all);

22 ab(2,all) -= (ab(2,0) / ab(0,0)) * ab(0,all);

23 ab(2,all) -= (ab(2,1) / ab(1,1)) * ab(1,all);

25 std::cout << "Matriz escalonada:\n" << ab << std::endl;

27 // redução

28 ab(1,all) -= (ab(1,2) / ab(2,2)) * ab(2,all);

29 ab(0,all) -= (ab(0,2) / ab(2,2)) * ab(2,all);

30 ab(0,all) -= (ab(0,1) / ab(1,1)) * ab(1,all);

31 ab(1,all) /= ab(1,1);

32 ab(2,all) /= ab(2,2);

33 ab(0,all) /= ab(0,0);

35 std::cout << "Matriz reduzida:\n" << ab << std::endl;

37 // solução

38 Eigen::VectorXd x(ab.rows());

39 x << ab(all, ab.cols()-1);

41 std::cout << "Solução x = \n" << x << std::endl;

42 return 0;

43}

4.4 Operações elemento-a-elemento

Código 19: ExOpElem.cpp

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 // mult. elemento-a-elemento

8 Eigen::ArrayXd arr {{1,2,3}};

9 Eigen::VectorXd vec {{4,5,6}};

11 std::cout << "arr * vec.array() =\n"

12 << arr * vec.array() << std::endl;

14 // divisão elemento-a-elemento

15 Eigen::MatrixXd a {{3, 6}, {12, 21}};

16 Eigen::MatrixXd b {{3, 2}, {2, 7}};

18 std::cout << "a.array() * b.array() =\n"

19 << a.array() / b.array() << std::endl;

21 // funções elementares

22 arr.resize(5);

23 arr << 0.0 , M_PI/3, M_PI_4,

24 M_PI/3, M_PI_2;

25 std::cout << "sin(arr) =\n"

26 << arr.sin() << std::endl;

28 std::cout << "b**2 =\n"

29 << b.array().pow(2) << std::endl;

31 return 0;

32}

arr * vec.array() =

a.array() * b.array() =

1 3

6 3

sin(arr) =

0.866025

0.707107

0.866025

b**2 =

9 4

4 49

Observação 4.4.1.(Mais sobre operações elemento-a-elemento)

Consulte mais sobre operações elemento-a-elemento em Eigen::The Array class and coefficient-wise operations.

Exercício 4.4.1.

Faça um código para computar os valores da função cosseno aplicada ao vetor

\boldsymbol{\theta}=(0,30^{\circ},45^{\circ},60^{\circ},90^{\circ}).

(24)

Resposta 0.

⬇

1#include <iostream>

2#include <cmath>

3#include <eigen3/Eigen/Eigen>

5int main()

7 Eigen::VectorXd theta {{0.0,

8 M_PI/6,

9 M_PI/4,

10 M_PI/3,

11 M_PI/2}};

12 std::cout << "cos(theta) = " << std::endl

13 << theta.array().cos() << std::endl;

15 return 0;

16}

Exercício 4.4.2.

Considere que os seguintes objetos estão declarados declarados:

⬇

1Eigen::ArrayXd arr {{1,2,3}};

2Eigen::VectorXd vec {{4,5,6}};

3Eigen::MatrixXd a {{3, 6, 4},

4 {2, 5, 3}};

Antes de implementar para verificar, forneça o resultado, quando definido, de cada uma das seguintes operações:

a)

arr - vec.array();
b)

a.array() * vec;
c)

a * vec;
d)

a.row(0).array() * arr;
e)

a.col(1) * vec;
f)

a.col(1).array() * arr;

Resposta 0.

a) $(-3,-3,-3)$ ; b) não definido; c) $(66,51)$ ; d) não definido.; e) $51$ ; f) não definido.

4.5 Inicialização

A biblioteca Eigen conta com métodos para inicialização de matrizes e arranjos.

•

Zero()

⬇

1arr = Eigen::ArrayXd::Zero(3)

⬇

arr =

0

0

0

⬇

1v = Eigen::RowVectorXd::Zero(2)

⬇

v = 0 0

⬇

1a = Eigen::MatrixXd::Zero(2,3)

⬇

a =

0 0 0

0 0 0
•

Ones()

⬇

1arr = Eigen::ArrayXd::Ones(2)

⬇

arr =

1

1
•

Constant()

⬇

1a = Eigen::MatrixXd::Constant(2,3,0.5)

⬇

a =

0.5 0.5 0.5

0.5 0.5 0.5
•

Random()

⬇

1a = Eigen::MatrixXd::Random(3,2)

⬇

a =

0.680375 0.59688

-0.211234 0.823295

0.566198 -0.604897
•

Identity()

⬇

1a = Eigen::MatrixXd::Identity(3,3)

⬇

a =

1 0 0

0 1 0

0 0 1
•

LinSpaced()

⬇

1arr = Eigen::RowVectorXd::LinSpaced(5, 0, 1)

⬇

arr =

0 0.25 0.5 0.75 1

Exercício 4.5.1.

Considere as seguintes partições uniformes de intervalos

	$\displaystyle x_{[0,1]}=\{0=x_{0}<x_{1}<\cdots<x_{n_{x}+1}=1\},$		(25)
	$\displaystyle y_{[1,2]}=\{1=y_{0}<y_{1}<\cdots<y_{n_{y}+1}=2\},$		(26)

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4void meshgrid(const Eigen::ArrayXd& x, const Eigen::ArrayXd& y,

5 Eigen::ArrayXXd& X, Eigen::ArrayXXd& Y)

7 X = x.replicate(1, y.size());

8 Y = y.transpose().replicate(x.size(), 1);

11int main()

12{

13 Eigen::ArrayXd x = Eigen::ArrayXd::LinSpaced(5, 0, 1);

14 Eigen::ArrayXd y = Eigen::ArrayXd::LinSpaced(3, 1, 2);

16 Eigen::ArrayXXd X;

17 Eigen::ArrayXXd Y;

18 meshgrid(x, y, X, Y);

20 std::cout << "X:\n" << X << "\n";

21 std::cout << "Y:\n" << Y << "\n";

23 return 0;

24}

Exercício 4.5.2.

Implemente uma função para computar a inverse de uma matriz não-singular pelo método de eliminação gaussiana. Verifique seu código com o método A.inverse() da Eigen.

Resposta 0.

⬇

1Eigen::MatrixXd inverse(const Eigen::MatrixXd& A)

3 // extendida

4 Eigen::MatrixXd E(A.rows(), 2 * A.cols());

5 E << A, Eigen::MatrixXd::Identity(A.rows(), A.cols());

7 // escalona

8 for (int j = 0; j < A.cols()-1; ++j) {

9 for (int i = j+1; i < A.rows(); ++i) {

10 E.row(i) -= E(i,j)/E(j,j) * E.row(j);

11 }

12 }

13 // reduzida

14 for (int j = A.cols()-1; j > 0; --j) {

15 for (int i = j-1; i >= 0; --i) {

16 E.row(i) -= E(i,j)/E(j,j) * E.row(j);

17 }

18 }

19 // normaliza

20 for (int j = 0; j < A.cols(); ++j) {

21 E.row(j) /= E(j,j);

22 }

24 return E.rightCols(A.cols());

25}

4.6 Reduções

Reduções são operações que transformam uma matriz ou um arranjo em um escalar.

•

sum()

⬇

1Eigen::MatrixXd a {{1,2,3},

2 {4,5,6}};

3double s = a.sum();

4std::cout << "s = " << s << std::endl;

⬇

s = 21
•

prod()

⬇

1double p = a.prod();

2std::cout << "p = " << p << std::endl;

⬇

p = 720
•

mean()

⬇

1double m = a.mean();

2std::cout << "m = " << m << std::endl;

⬇

m = 3.5
•

minCoeff()/maxCoeff()

⬇

1double min = a.minCoeff();

2std::cout << "min = " << min << std::endl;

3double max = a.maxCoeff();

4std::cout << "max = " << max << std::endl;

⬇

min = 1

max = 6
•

trace()

⬇

1Eigen::MatrixXd b {{1,3},

2 {4,6}};

3

4double tr = b.trace();

5std::cout << "tr = " << tr << std::endl;

⬇

tr = 7
•

norm()

⬇

1Eigen::VectorXd v {{1,2,3}};

2Eigen::MatrixXd m {{1,2},

3 {3,4}};

4// norma l2

5double v_norm = v.norm();

6std::cout << "v_norm = " << v_norm << std::endl;

7// norma L2 (Frobenius)

8double m_norm = m.norm();

9std::cout << "m_norm = " << m_norm << std::endl;

⬇

n = 3.74166

m_norm = 5.47723

Exercício 4.6.1.

Faça um código que aloque os seguintes vetores

	$\displaystyle\vec{u}=(1,2,3),$		(28)
	$\displaystyle\vec{v}=(4,5,6),$		(29)

e compute a o produto interno (ou escalar) $u\cdot v$ . Verifique seu resultado com o método u.dot(v) da Eigen.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::VectorXd u {{1,2,3}};

7 Eigen::VectorXd v {{4,5,6}};

9 Eigen::VectorXd w(u.size());

10 w << u.array() * v.array();

11 std::cout << "u.v = " << w.sum() << std::endl;

12 std::cout << "u.dot(v) = " << u.dot(v) << std::endl;

14 return 0;

15}

Exercício 4.6.2.

Implemente uma função para computar o erro médio quadrático entre dois vetores u e v de mesmo tamanho. O erro médio quadrático é dado por

\text{EMQ}=\frac{1}{n}\sum_{i=1}^{n}(u_{i}-v_{i})^{2}.

(30)

Resposta 0.

⬇

1double EMQ(const Eigen::VectorXd& u, const Eigen::VectorXd& v)

3 return ((u - v).array().pow(2)).sum() / u.size();

4.7 Sistemas lineares

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

usando a o método da decomposição LU com pivotamento parcial.

Código 20: ExSistLin.cpp

⬇

1Eigen::MatrixXd m {{2,-1},

2 {1,3}};

3Eigen::VectorXd b {{3,-2}};

5std::cout << "x =\n"

6 << m.lu().solve(b) << std::endl;

x =

-1

Observação 4.7.1.(Mais sobre sistemas lineares)

Consulte mais sobre sistemas lineares em Eigen::Linear Algebra.

Exercício 4.7.1.

Considere o seguinte sistema de equações lineares

$\displaystyle 2x_{1}+3x_{2}+4x_{3}$	$\displaystyle=-1,$	(33)
$\displaystyle 3x_{1}-4x_{2}+5x_{3}$	$\displaystyle=7,$	(34)
$\displaystyle x_{1}+x_{2}-x_{3}$	$\displaystyle=0.$	(35)

Escreva um código que aloque a matriz dos coeficientes $A$ e o vetor dos termos constantes $\vec{b}$ do sistema. Então, compute a solução do sistema usando o método LU com pivotamento parcial.

Resposta 0.

⬇

1#include <iostream>

2#include <eigen3/Eigen/Eigen>

4int main()

6 Eigen::MatrixXd a {{2,3,4},

7 {3,-4,5},

8 {1,1,-1}};

9 Eigen::VectorXd b {{-1,7,0}};

11 Eigen::VectorXd x = a.lu().solve(b);

12 std::cout << "x =\n" << x << std::endl;

14 return 0;

15}

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

| | | |