Bastão de Asclépio & Símbolo de Integral
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL","pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(expm, warn.conflicts=FALSE))
suppressMessages(library(far, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(igraph, warn.conflicts=FALSE))
suppressMessages(library(igraphdata, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(matlib, warn.conflicts=FALSE))
suppressMessages(library(Matrix, warn.conflicts=FALSE))
suppressMessages(library(matrixcalc, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(rayshader, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
Desmos Matrix: https://www.desmos.com/matrix?lang=pt-BR
Desmos 2D Modo Complexo: https://www.desmos.com/calculator?lang=pt-BR
Matrix Approach to Optimization: Hessian Determinant: Dr Bilal Mehmood
https://rich-iannone.github.io/DiagrammeR/
DiagrammeR::DiagrammeR
Algebra linear aplicada em Análise de Teste Diagnóstico:
Exemplo 14.3.4: Genética: Exploration 4.4: A Genetics Model, p. 289:
Linear Algebra and Its Applications with R - Yoshida - 2021.pdf
magick
package in R to demonstrate image processing
without getting into detailslpSolve
packagegMOIP
packageigraph
packageigraphdata
package4.4. Matrix Models of Base Substitution, p 138: “We begin by modeling the ancestral sequence probabilistically. Each site in the sequence is one of the four bases A, G, C, or T , chosen randomly according to some probabilities”:
Basics of Matrix Algebra for Statistics with R - Fieller - 2015.pdf
Matrix Operations - Bronson - 1989.pdf
RPubs
10. Matriz, vetor e número complexo (Capítulos 14 e 15)
A teoria de matrizes estuda as propriedades, operações e aplicações de matrizes, que são arranjos retangulares de números dispostos em linhas e colunas. Formalmente, uma matriz \(A \in \mathbb{R}^{m \times n}\) é um objeto com \(m\) linhas e \(n\) colunas:
\[ A = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix} \]
Operações básicas:
Aplicações:
A teoria de matrizes é a base da álgebra linear, com papel fundamental na matemática aplicada, estatística, engenharia e ciência de dados.
Matrix (mathematics): Wikipedia
Nas últimas décadas as matrizes tornaram-se indispensáveis em diversas aplicações da matemática à biologia, principalmente na estatística.
Um caso especial importante de matrizes são os vetores. Vamos nos concentrar em vetores e suas diversas aplicações nas ciências da vida.
Matriz é um arranjo retangular de objetos do mesmo tipo.
[,1] [,2] [,3]
[1,] 2 1 -1
[2,] -3 -1 2
[3,] -2 1 2
[1] 2 3 -1
# provide implicit or explicit labels
matlib::printMatEqn("Matriz A" = A, "*", "Vetor x" = x,
"=", "Vetor b" = A %*% x)
Matriz A Vetor x Vetor b
2 1 -1 * 2 = 8
-3 -1 2 3 -11
-2 1 2 -1 -3
A x b
2 1 -1 * 2 = 8
-3 -1 2 3 -11
-2 1 2 -1 -3
A x A %*% x
2 1 -1 * 2 = 8
-3 -1 2 3 -11
-2 1 2 -1 -3
A x b
2 1 -1 x1 = 4
-3 -1 2 x2 2
-2 1 2 x3 1
2*x1 + 1*x2 - 1*x3 = 4
-3*x1 - 1*x2 + 2*x3 = 2
-2*x1 + 1*x2 + 2*x3 = 1
# decimal example
A <- matrix(c(0.5, 1, 3, 0.75, 2.8, 4),
nrow=2, ncol=3, byrow=TRUE)
x <- c(0.5, 3.7, 2.3)
y <- c(0.7, -1.2)
b <- A %*% x - y
b
[,1]
[1,] 10.150
[2,] 21.135
A x y b
0.50 1.00 3.00 * 0.5 - 0.7 = 10.150
0.75 2.80 4.00 3.7 -1.2 21.135
2.3
A x y b
1/2 1 3 * 1/2 - 7/10 = 203/20
3/4 14/5 4 37/10 -6/5 4227/200
23/10
As matrizes são escritas entre parênteses ou colchetes.
Com barras verticais não representa uma matriz, mas a quantidade \(ad - bc\). Este símbolo é chamado de determinante.
\[ \det\left[\begin{array}{rr} a & b\\ c & d \end{array}\right]=\left|\begin{array}{rr} a & b\\ c & d \end{array}\right|=ad - bc \]
Uma notação útil são os índices duplos. Suponha que um pesquisador realize um experimento repetidamente sob várias condições usando diferentes tratamentos. Denotamos o número de tratamentos por \(m\) e o número de medições em cada tratamento por \(n\). Se \(n\) e \(m\) não forem particularmente pequenos, o alfabeto não conteria letras suficientes para denotar as medidas individuais, nem subscritos simples como em \(x_1,x_2,\ldots\) seriam úteis. É conveniente introduzir um segundo subscrito.
Assim, por \(x_{23}\) (leia-se: \(x\) dois, três) queremos dizer a terceira medição do segundo tratamento. As medidas formam a matriz
\[ \begin{array}{c c c c c c c} & {\text{Medida}} \\ & & 1 & 2 & 3 & \cdots & n \\ \text{Tratamento} & 1 & x_{11} & x_{12} & x_{13} & \cdots & x_{1n} \\ & 2 & x_{21} & x_{22} & x_{23} & \cdots & x_{2n} \\ & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ & m & x_{m1} & x_{m2} & x_{m3} & \cdots & x_{mn} \\ \end{array} \]
Cada \(x\) é chamado de elemento ou componente da matriz. O elemento geral na \(i\)-ésima linha e nas \(j\)-ésimas colunas é \(x_{ij}\). O primeiro subscrito designa a linha, o segundo a coluna na qual o elemento está localizado. Em nossa notação, \(i\) varia de 1 a \(m\), e \(j\) varia de 1 a \(n\).
Por brevidade, as matrizes são denotadas por letras únicas. Isso pode causar algum mal-entendido, pois as letras frequentemente representam números únicos e não matrizes de números. Portanto, alguma distinção é desejável.
Na livro adotado, são usadas letras em negrito como notações para matrizes. Por exemplo,
\[ \textbf{A}=\left[\begin{array}{rrr} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23} \end{array}\right] \]
Para maior clareza, o número de linhas e colunas geralmente é adicionado à notação. Assim, \(A\) é uma matriz “2 por 3”, simbolicamente \(\underset{2\times3}{\text{A}}\).
Neste texto, não adotaremos negrito como notação para matriz por ser redundante com o número de linhas e colunas.
Um vetor é uma matriz com o número de linhas ou colunas igual a 1.
Então \(u_{1\times 4}\) é um vetor linha e \(\underset{4\times 1}{u}\) é um vetor coluna. Note que um vetor linha transposto é um vetor coluna, i.e., \(\underset{1\times 4}{u^{\prime}}=\underset{4\times 1}{u}\)
Transpose[{{a, b, c}, {d, e, f}}]
[,1] [,2] [,3]
[1,] "a" "b" "c"
[,1] [,2] [,3]
[1,] a b c
[,1] [,2] [,3]
[1,] "a" "c" "e"
[2,] "b" "d" "f"
[,1] [,2] [,3]
[1,] a c e
[2,] b d f
[1] 1 3
[1] "matrix" "array"
[1] 2 3
[1] "matrix" "array"
[,1]
[1,] "a"
[2,] "b"
[3,] "c"
[1] 3 1
[,1] [,2]
[1,] "a" "b"
[2,] "c" "d"
[3,] "e" "f"
[1] 3 2
Um sistema linear consiste em um conjunto de equações lineares com as mesmas variáveis. A forma geral de um sistema linear com \(m\) equações e \(n\) incógnitas é:
\[ \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\ \quad\vdots \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n = b_m \end{cases} \]
ou, na forma matricial:
\[ \mathbf{A}\mathbf{x} = \mathbf{b} \]
onde: - \(A \in \mathbb{R}^{m \times n}\) é a matriz dos coeficientes, - \(\mathbf{x} \in \mathbb{R}^n\) é o vetor incógnita, - \(\mathbf{b} \in \mathbb{R}^m\) é o vetor dos termos independentes.
Um sistema linear pode ser classificado quanto à existência e unicidade das soluções:
Existe uma única solução. O sistema é compatível e tem solução única. Isso ocorre quando:
\[ \text{posto}(A) = \text{posto}(A|b) = n \]
Existem infinitas soluções. Ocorre quando:
\[ \text{posto}(A) = \text{posto}(A|b) < n \]
Não há solução. Isso ocorre quando:
\[ \text{posto}(A) < \text{posto}(A|b) \]
Essas classificações são baseadas no Teorema de Rouché-Capelli e podem ser verificadas via cálculo do posto (rank) das matrizes envolvidas.
# Usar função `qr` para verificar posto
qr_rank <- function(M) qr(M)$rank
# Exemplo 1: Sistema Possível e Determinado (SPD)
A1 <- matrix(c(2, 1,
1, 1), nrow = 2, byrow = TRUE)
b1 <- c(4, 3)
solve(A1, b1) # solução única
[1] 1 2
rA1 <- qr_rank(A1)
rA1b <- qr_rank(cbind(A1, b1))
c(rA1, rA1b) # mesmo posto e igual ao número de variáveis
[1] 2 2
# Exemplo 2: Sistema Possível e Indeterminado (SPI)
A2 <- matrix(c(1, 2,
2, 4), nrow = 2, byrow = TRUE)
b2 <- c(6, 12)
try(solve(A2, b2))
Error in solve.default(A2, b2) :
Lapack routine dgesv: system is exactly singular: U[2,2] = 0
rA2 <- qr_rank(A2)
rA2b <- qr_rank(cbind(A2, b2))
c(rA2, rA2b) # mesmo posto, mas menor que número de variáveis
[1] 1 1
# Exemplo 3: Sistema Impossível (SI)
A3 <- matrix(c(1, 2,
2, 4), nrow = 2, byrow = TRUE)
b3 <- c(6, 13)
try(solve(A3, b3))
Error in solve.default(A3, b3) :
Lapack routine dgesv: system is exactly singular: U[2,2] = 0
rA3 <- qr_rank(A3)
rA3b <- qr_rank(cbind(A3, b3))
c(rA3, rA3b) # postos diferentes, sistema é impossível
[1] 1 2
O sistema linear a seguir é consistente e determinado. Portanto, tem solução única.
\[ \begin{align} x-5y&=2\\ 3x+4y&=8 \end{align} \]
\[ \left[\begin{array}{rrr} 1 & -5\\ 3 & 4 \end{array}\right] \left[\begin{array}{r} x\\ y \end{array}\right]= \left[\begin{array}{r} 2\\ 8 \end{array}\right] \]
A primeira matriz contém os coeficientes conhecidos de \(x\) e \(y\) em sua ordem apropriada. A segunda matriz contém as incógnitas e a última matriz as quantidades conhecidas no lado direito das equações.
A solução é \(x=\dfrac{48}{19}\approx 2.526\) e \(y=\dfrac{2}{19}\approx 0.105\).
solve x - 5 y = 2, 3 x + 4 y = 8
Solve[{x - 5 y == 2, 3 x + 4 y == 8}, {x, y}]
LinearSolve[{{1, -5}, {3, 4}}, {2, 8}]
[,1] [,2]
[1,] 1 -5
[2,] 3 4
[1] 2 8
1*x1 - 5*x2 = 2
3*x1 + 4*x2 = 8
x[1] - 5*x[2] = 2
3*x[1] + 4*x[2] = 8
[1] 2.5263158 0.1052632
[1] 48/19 2/19
[,1]
[1,] 2.5263158
[2,] 0.1052632
[,1]
[1,] 2.5263158
[2,] 0.1052632
1*48/19 - 5*2/19 = 2
3*48/19 + 4*2/19 = 8
x1 - 5*x2 = 2
3*x1 + 4*x2 = 8
\begin{array}{lllll}
1 \cdot x_1 &-& 5 \cdot x_2 &=& 2 \\
3 \cdot x_1 &+& 4 \cdot x_2 &=& 8 \\
\end{array}
a=[1 -5;3 4];
b=[2;8];
linsolve(a,-b)
//Método direto de solução do sistema linear por meio da equação ax-b=0.
//Notar que b deve conter apenas valores não negativos, i.e., b>=0.
quit
Scilab 2023.1.0 (May 23 2023, 09:23:00)
ans =
2.5263158
0.1052632
O sistema linear a seguir é consistente e determinado. Portanto, tem solução única.
\[ \begin{align} 13x-4y+2z&=1\\ 4x+11y-2z&=2\\ 2x-2y+8z&=4 \end{align} \]
{{13, -4, 2}, {4, 11, -2}, {2, -2, 8}} . {x, y, z} = {1, 2, 4}
# Sistema linear determinado consistente: 3D
A <- matrix(c(13, -4, 2,
4, 11, -2,
2, -2, 8),
byrow=TRUE, nrow=3)
A
[,1] [,2] [,3]
[1,] 13 -4 2
[2,] 4 11 -2
[3,] 2 -2 8
13*x1 - 4*x2 + 2*x3 = 1
4*x1 + 11*x2 - 2*x3 = 2
2*x1 - 2*x2 + 8*x3 = 4
[1] 0.07142857 0.25510204 0.54591837
[1] 1/14 25/98 107/196
13*1/14 - 4*25/98 + 2*107/196 = 1
4*1/14 + 11*25/98 - 2*107/196 = 2
2*1/14 - 2*25/98 + 8*107/196 = 4
13*x1 - 4*x2 + 2*x3 = 1
4*x1 + 11*x2 - 2*x3 = 2
2*x1 - 2*x2 + 8*x3 = 4
\begin{array}{lllllll}
13 \cdot x_1 &-& 4 \cdot x_2 &+& 2 \cdot x_3 &=& 1 \\
4 \cdot x_1 &+& 11 \cdot x_2 &-& 2 \cdot x_3 &=& 2 \\
2 \cdot x_1 &-& 2 \cdot x_2 &+& 8 \cdot x_3 &=& 4 \\
\end{array}
# Criar a matriz aumentada [A|b]
Ab <- cbind(A, b)
# Obter a forma escalonada por linhas
echelon_result <- matlib::echelon(Ab, verbose = TRUE)
Initial matrix:
b
[1,] 13 -4 2 1
[2,] 4 11 -2 2
[3,] 2 -2 8 4
row: 1
multiply row 1 by 0.07692308
b
[1,] 1 -0.3076923 0.1538461 0.07692308
[2,] 4 11.0000000 -2.0000000 2.00000000
[3,] 2 -2.0000000 8.0000000 4.00000000
multiply row 1 by 4 and subtract from row 2
b
[1,] 1 -0.3076923 0.1538461 0.07692308
[2,] 0 12.2307692 -2.6153846 1.69230769
[3,] 2 -2.0000000 8.0000000 4.00000000
multiply row 1 by 2 and subtract from row 3
b
[1,] 1 -0.3076923 0.1538461 0.07692308
[2,] 0 12.2307692 -2.6153846 1.69230769
[3,] 0 -1.3846154 7.6923077 3.84615385
row: 2
multiply row 2 by 0.08176101
b
[1,] 1 -0.3076923 0.1538461 0.07692308
[2,] 0 1.0000000 -0.2138365 0.13836478
[3,] 0 -1.3846154 7.6923077 3.84615385
multiply row 2 by 0.3076923 and add to row 1
b
[1,] 1 0.000000 0.08805031 0.1194969
[2,] 0 1.000000 -0.21383648 0.1383648
[3,] 0 -1.384615 7.69230769 3.8461538
multiply row 2 by 1.384615 and add to row 3
b
[1,] 1 0 0.08805031 0.1194969
[2,] 0 1 -0.21383648 0.1383648
[3,] 0 0 7.39622642 4.0377358
row: 3
multiply row 3 by 0.1352041
b
[1,] 1 0 0.08805031 0.1194969
[2,] 0 1 -0.21383648 0.1383648
[3,] 0 0 1.00000000 0.5459184
multiply row 3 by 0.08805031 and subtract from row 1
b
[1,] 1 0 0.0000000 0.07142857
[2,] 0 1 -0.2138365 0.13836478
[3,] 0 0 1.0000000 0.54591837
multiply row 3 by 0.2138365 and add to row 2
b
[1,] 1 0 0 0.07142857
[2,] 0 1 0 0.25510204
[3,] 0 0 1 0.54591837
# Verificar a inconsistência
# Verificamos se há uma linha onde todos os coeficientes são zero, mas o termo constante não é zero
is_inconsistent <- any(rowSums(echelon_result[, 1:ncol(A)]) == 0 &
echelon_result[, ncol(A) + 1] != 0)
if (is_inconsistent) {
cat("O sistema é inconsistente.\n")
} else {
cat("O sistema é consistente.\n")
}
O sistema é consistente.
Numa fábrica de uma multinacional farmacêutica os últimos estágios na produção de xaropes são a esterilização das embalagens, o enchimento, a etiquetagem e o embalamento.
A empresa fabrica quatro xaropes diferentes: A, B, C e D.
O tempo gasto por embalagem em cada etapa é:
Xarope | Esterilização (min) | Enchimento (min) | Etiquetagem (min) | Embalamento (min) |
---|---|---|---|---|
A | 2 | 0.5 | 0.25 | 3 |
B | 2 | 0.5 | 0.25 | 4 |
C | 2.25 | 0.75 | 0.25 | 5 |
D | 1.75 | 0.5 | 0.25 | 3 |
As máquinas disponíveis têm os seguintes limites diários de tempo:
Máquina | Tempo máximo (min) |
---|---|
Esterilização | 520 |
Enchimento | 146 |
Etiquetagem | 65.75 |
Embalamento | 955 |
Pergunta:
Quantas embalagens de cada tipo de xarope devem ser produzidas por dia, respeitando os tempos máximos das máquinas?
Modelagem:
Seja:
As restrições são dadas pelas seguintes equações:
\[ \begin{cases} 2x_1 + 2x_2 + 2.25x_3 + 1.75x_4 = 520 \quad \text{(Esterilização)} \\ 0.5x_1 + 0.5x_2 + 0.75x_3 + 0.5x_4 = 146 \quad \text{(Enchimento)} \\ 0.25x_1 + 0.25x_2 + 0.25x_3 + 0.25x_4 = 65.75 \quad \text{(Etiquetagem)} \\ 3x_1 + 4x_2 + 5x_3 + 3x_4 = 955 \quad \text{(Embalamento)} \end{cases} \]
A <- matrix(c(
2, 2, 2.25, 1.75,
0.5, 0.5, 0.75, 0.5,
0.25, 0.25, 0.25, 0.25,
3, 4, 5, 3
), nrow = 4, byrow = TRUE)
b <- c(520, 146, 65.75, 955)
solve(A, b)
[1] 73 50 58 82
A solução é \((x_1,x_2,x_3,x_4)=(73, 50, 58, 82)\).
Assim, para duas equações lineares simultâneas com três incógnitas (sistema subdeterminado) como
\[ \begin{align} a_1x+b_1y+c_1z&=d_1\\ a_2x+b_2y+c_2z&=d_2 \end{align} \]
podemos representar o sistema com três matrizes
\[ \left[\begin{array}{rrr} a_1 & b_1 & c_1\\ a_2 & b_2 & c_2 \end{array}\right] \left[\begin{array}{r} x\\ y\\ z \end{array}\right]= \left[\begin{array}{r} d_1\\ d_2 \end{array}\right] \]
O sistema
\[ \begin{align} 2x+2y-z&=1\\ 2x-y+z&=2 \end{align} \]
tem solução
\[ \left\{(x,y,z):x=\dfrac{5}{6}-\dfrac{1}{6}z, \; y=-\dfrac{1}{3}+\dfrac{2}{3}z\right\} \]
ou
\[ \left\{(x,y,z):y=3-4x, \; z=5-6x\right\} \]
{{a_1, b_1, c_1}, {a_2, b_2, c_2}} . {x, y, z} = {d_1, d_2}
{{2, 2, -1}, {2, -1, 1}} . {x, y, z} = {1, 2}
solve 2x + 2y - z = 1, 2x - y + z = 2
LinearSolve
resolve de maneira incompletaSistema sobredeterminado ocorre quando há mais equações que incógnitas ($m > n$). Isso pode levar a duas situações:
Sistema sobredeterminado consistente: todas as equações se interceptam num mesmo ponto, ou seja, há pelo menos uma solução que satisfaz todas simultaneamente.
Sistema sobredeterminado inconsistente: não há solução que satisfaça todas as equações ao mesmo tempo — por exemplo, quando as equações representam retas (ou planos) que não se encontram num único ponto comum.
Exemplo inconsistente:
\[ \begin{cases} 2x - y = 1 \\ x + y = 2 \\ x - y = 3 \end{cases} \]
Esse sistema possui três equações e duas incógnitas. Não há solução $(x, y)$ que satisfaça as três ao mesmo tempo. A matriz aumentada na forma escalonada gera uma linha como:
\[ 0x + 0y = c \quad \text{com } c \ne 0 \]
o que indica contradição e, portanto, inconsistência.
O teste programado no R com matlib::echelon()
detecta
essa inconsistência.
Quando o sistema é inconsistente, pode-se usar mínimos quadrados:
\[ \widehat{y} = (A^{\prime} A)^{-1} A^{\prime} b \]
que fornece a melhor aproximação possível da solução, minimizando a soma dos quadrados dos resíduos.
O sistema linear sobredeterminado a seguir é consistente e, portanto, tem solução.
\[ \begin{align} 2x+2y&=1\\ x-y&=2\\ 3x+y&=3 \end{align} \]
{{2, 2}, {1, -1}, {3, 1}} . {x, y} = {1, 2, 3}
# Sistema linear sobredeterminado consistente
A <- matrix(c(2, 2,
1, -1,
3, 1),
byrow=TRUE, nrow=3)
A
[,1] [,2]
[1,] 2 2
[2,] 1 -1
[3,] 3 1
2*x1 + 2*x2 = 1
1*x1 - 1*x2 = 2
3*x1 + 1*x2 = 3
2*x[1] + 2*x[2] = 1
x[1] - 1*x[2] = 2
3*x[1] + x[2] = 3
Error in solve.default(A, b) : 'a' (3 x 2) deve ser quadrada
[,1]
[1,] 1.25
[2,] -0.75
[,1]
[1,] 5/4
[2,] -3/4
2*5/4 + 2*-3/4 = 1
1*5/4 - 1*-3/4 = 2
3*5/4 + 1*-3/4 = 3
2*x1 + 2*x2 = 1
x1 - 1*x2 = 2
3*x1 + x2 = 3
\begin{array}{lllll}
2 \cdot x_1 &+& 2 \cdot x_2 &=& 1 \\
1 \cdot x_1 &-& 1 \cdot x_2 &=& 2 \\
3 \cdot x_1 &+& 1 \cdot x_2 &=& 3 \\
\end{array}
# Criar a matriz aumentada [A|b]
Ab <- cbind(A, b)
# Obter a forma escalonada por linhas
echelon_result <- matlib::echelon(Ab, verbose = TRUE)
Initial matrix:
b
[1,] 2 2 1
[2,] 1 -1 2
[3,] 3 1 3
row: 1
exchange rows 1 and 3
b
[1,] 3 1 3
[2,] 1 -1 2
[3,] 2 2 1
multiply row 1 by 0.3333333
b
[1,] 1 0.3333333 1
[2,] 1 -1.0000000 2
[3,] 2 2.0000000 1
subtract row 1 from row 2
b
[1,] 1 0.3333333 1
[2,] 0 -1.3333333 1
[3,] 2 2.0000000 1
multiply row 1 by 2 and subtract from row 3
b
[1,] 1 0.3333333 1
[2,] 0 -1.3333333 1
[3,] 0 1.3333333 -1
row: 2
exchange rows 2 and 3
b
[1,] 1 0.3333333 1
[2,] 0 1.3333333 -1
[3,] 0 -1.3333333 1
multiply row 2 by 0.75
b
[1,] 1 0.3333333 1.00
[2,] 0 1.0000000 -0.75
[3,] 0 -1.3333333 1.00
multiply row 2 by 0.3333333 and subtract from row 1
b
[1,] 1 0.000000 1.25
[2,] 0 1.000000 -0.75
[3,] 0 -1.333333 1.00
multiply row 2 by 1.333333 and add to row 3
b
[1,] 1 0 1.25
[2,] 0 1 -0.75
[3,] 0 0 0.00
row: 3
# Verificar a inconsistência
# Verificamos se há uma linha onde todos os coeficientes são zero, mas o termo constante não é zero
is_inconsistent <- any(rowSums(echelon_result[, 1:ncol(A)]) == 0 &
echelon_result[, ncol(A) + 1] != 0)
if (is_inconsistent) {
cat("O sistema é inconsistente.\n")
} else {
cat("O sistema é consistente.\n")
}
O sistema é consistente.
O sistema linear sobredeterminado a seguir é inconsistente e, portanto, não tem solução.
\[ \begin{align} 2x-y&=1\\ x+y&=2\\ x-y&=3 \end{align} \]
{{2, -1}, {1, 1}, {1, -1}} . {x, y} = {1, 2, 3}
# Sistema linear sobredeterminado inconsistente
A <- matrix(c(2, -1,
1, 1,
1, -1),
byrow=TRUE, nrow=3)
A
[,1] [,2]
[1,] 2 -1
[2,] 1 1
[3,] 1 -1
2*x1 - 1*x2 = 1
1*x1 + 1*x2 = 2
1*x1 - 1*x2 = 3
2*x[1] - 1*x[2] = 1
x[1] + x[2] = 2
x[1] - 1*x[2] = 3
Error in solve.default(A, b) : 'a' (3 x 2) deve ser quadrada
[,1]
[1,] 1.2142857
[2,] 0.1428571
[,1]
[1,] 17/14
[2,] 1/7
2*17/14 - 1*1/7 = 1
1*17/14 + 1*1/7 = 2
1*17/14 - 1*1/7 = 3
2*x1 - 1*x2 = 1
x1 + x2 = 2
x1 - 1*x2 = 3
\begin{array}{lllll}
2 \cdot x_1 &-& 1 \cdot x_2 &=& 1 \\
1 \cdot x_1 &+& 1 \cdot x_2 &=& 2 \\
1 \cdot x_1 &-& 1 \cdot x_2 &=& 3 \\
\end{array}
# Criar a matriz aumentada [A|b]
Ab <- cbind(A, b)
# Obter a forma escalonada por linhas
echelon_result <- matlib::echelon(Ab, verbose = TRUE)
Initial matrix:
b
[1,] 2 -1 1
[2,] 1 1 2
[3,] 1 -1 3
row: 1
multiply row 1 by 0.5
b
[1,] 1 -0.5 0.5
[2,] 1 1.0 2.0
[3,] 1 -1.0 3.0
subtract row 1 from row 2
b
[1,] 1 -0.5 0.5
[2,] 0 1.5 1.5
[3,] 1 -1.0 3.0
subtract row 1 from row 3
b
[1,] 1 -0.5 0.5
[2,] 0 1.5 1.5
[3,] 0 -0.5 2.5
row: 2
multiply row 2 by 0.6666667
b
[1,] 1 -0.5 0.5
[2,] 0 1.0 1.0
[3,] 0 -0.5 2.5
multiply row 2 by 0.5 and add to row 1
b
[1,] 1 0.0 1.0
[2,] 0 1.0 1.0
[3,] 0 -0.5 2.5
multiply row 2 by 0.5 and add to row 3
b
[1,] 1 0 1
[2,] 0 1 1
[3,] 0 0 3
row: 3
multiply row 3 by 0.3333333
b
[1,] 1 0 1
[2,] 0 1 1
[3,] 0 0 1
subtract row 3 from row 1
b
[1,] 1 0 0
[2,] 0 1 1
[3,] 0 0 1
subtract row 3 from row 2
b
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
# Verificar a inconsistência
# Verificamos se há uma linha onde todos os coeficientes são zero, mas o termo constante não é zero
is_inconsistent <- any(rowSums(echelon_result[, 1:ncol(A)]) == 0 &
echelon_result[, ncol(A) + 1] != 0)
if (is_inconsistent) {
cat("O sistema é inconsistente.\n")
} else {
cat("O sistema é consistente.\n")
}
O sistema é inconsistente.
Um sistema linear sobredeterminado consistente homogêneo é da forma:
\[ A\mathbf{x} = \mathbf{0} \]
com $A$ sendo uma matriz \(m \times n\) tal que \(m > n\), ou seja, há mais equações do que incógnitas. Como o sistema é homogêneo, ele sempre admite a solução trivial \(\mathbf{x} = \mathbf{0}\).
Se o posto de \(A\) for menor que \(n\), o sistema tem infinitas soluções. Caso o posto de \(A\) seja igual a \(n\), a única solução é a trivial.
Exemplo:
\[ \begin{cases} x + y = 0 \\ 2x + 2y = 0 \\ 3x + 3y = 0 \end{cases} \]
Esse sistema tem 3 equações e 2 incógnitas. É sobredeterminado, homogêneo e consistente, pois todas as equações são linearmente dependentes. A solução geral é \(x = -y\).
solve {x + y = 0, 2x + 2y = 0, 3x + 3y = 0}
[,1]
[1,] 1
[2,] -1
Um sistema linear homogêneo subdeterminado (menos equações que variáveis) é sempre:
Como geralmente ocorre com sistemas homogêneos de balanceamento químico, em que há mais substâncias do que equações de conservação, o sistema é:
Uma das áreas científicas com maior peso nas ciências farmacêuticas é a química.
A química estuda as transformações da matéria através de reações químicas. Para interpretar uma reação química é preciso balanceá-la.
O princípio subjacente a este balanceamento é o da conservação do número de átomos de cada elemento que diz que, para cada elemento, o número de átomos que existe nos reagentes tem que ser igual ao que aparece nos produtos.
A Lei de Lavoisier, ou Lei da Conservação da Massa, afirma que:
“Na natureza, nada se cria, nada se perde, tudo se transforma.”
Ou, em termos químicos:
“Em uma reação química realizada em sistema fechado, a massa total dos reagentes é igual à massa total dos produtos.”
Essa lei foi formulada por Antoine Laurent de Lavoisier no século XVIII e é a base do balanceamento das equações químicas.
Consequência prática:
O número de átomos de cada elemento deve ser conservado antes e depois da reação.
Qualquer estudante que tenha frequentado química no ensino pré-universitário não tem qualquer dificuldade em acertar mentalmente a equação
\[ \text{H}_2 + \text{O}_2 \rightarrow \text{H}_2\text{O} \]
colocando um 2 atrás da molécula de hidrogênio e da molécula de água:
\[ 2\,\text{H}_2 + \text{O}_2 \rightarrow 2\,\text{H}_2\text{O} \]
Considere a reação de combustão do etanol (\(\mathrm{C_2H_5OH}\)):
\[ \mathrm{C_2H_5OH} + \mathrm{O_2} \rightarrow \mathrm{CO_2} + \mathrm{H_2O} \]
Para balancear, igualamos o número de átomos de cada elemento nos dois lados.
Montamos o sistema linear a partir da conservação dos átomos:
Equação balanceada:
\[ a\,\mathrm{C_2H_5OH} + b\, \mathrm{O_2} \rightarrow c\, \mathrm{CO_2} + d\, \mathrm{H_2O} \]
Resolvemos o sistema homogêneo \(Ax = 0\) com \(x = (a, b, c, d)\):
solve 2a - c = 0, 6a - 2d = 0, a + 2b - 2c - d = 0 where a = 1
M <- matrix(c(
2, 0, -1, 0, # 2a - c = 0
6, 0, 0, -2, # 6a - 2d = 0
1, 2, -2, -1 # a + 2b - 2c - d = 0
), nrow=3, byrow=TRUE)
N <- pracma::nullspace(M)
round(N / min(N))
[,1]
[1,] 1
[2,] 3
[3,] 2
[4,] 3
A solução é \(x = (a, b, c, d)=(1,3,2,3)\).
Equação balanceada:
\[ \mathrm{C_2H_5OH} + 3\, \mathrm{O_2} \rightarrow 2\, \mathrm{CO_2} + 3\, \mathrm{H_2O} \]
A reação
\[ \text{As}_2\text{S}_5 + \text{HNO}_3 \rightarrow \text{H}_3\text{AsO}_4 + \text{H}_2\text{SO}_4 + \text{H}_2\text{O} + \text{NO}_2 \]
é uma reação de oxirredução entre o dissulfeto de diarsênio (As₂S₅) e o ácido nítrico (HNO₃), que atua como agente oxidante.
Pode ser classificada como:
– reação de oxidação de arseneto e sulfeto por ácido nítrico concentrado – reação redox inorgânica com liberação de NO₂ – reação típica de decomposição oxidativa com formação de oxoácidos (H₃AsO₄ e H₂SO₄)
Não possui nome tradicional específico na literatura, sendo corretamente descrita como oxidação de As₂S₅ por HNO₃ concentrado.
O princípio utilizado no balanceamento é sempre o mesmo – o número de
átomos de cada elemento conserva-se – porém, o balanceamento por
inspeção visual não é simples neste caso e a saída mais fácil é mesmo
recorrer à matemática, sendo que aquilo que procuramos é o valor dos
chamados coeficientes estequiométricos.
São as nossas incógnitas que representaremos pelas letras \(a\), \(b\), \(c\), \(d\), \(e\)
e \(f\):
\[ a\,\text{As}_2\text{S}_5 + b\,\text{HNO}_3 \rightarrow c\,\text{H}_3\text{AsO}_4 + d\,\text{H}_2\text{SO}_4 + e\,\text{H}_2\text{O} + f\,\text{NO}_2 \]
Montamos as equações de conservação para cada elemento:
solve 2a = c, 5a = d, b = 3c + 2d + 2e, b = f, 3b = 4c + 4d + e + 2f where a = 1
Reorganizamos as equações para a forma \(M x = 0\), com \(x = (a, b, c, d, e, f)^{\prime}\):
\[ \begin{bmatrix} 2 & 0 & -1 & 0 & 0 & 0 \\ 5 & 0 & 0 & -1 & 0 & 0 \\ 0 & 1 & -3 & -2 & -2 & 0 \\ 0 & 1 & 0 & 0 & 0 & -1 \\ 0 & 3 & -4 & -4 & -1 & -2 \\ \end{bmatrix} \cdot \begin{bmatrix} a \\ b \\ c \\ d \\ e \\ f \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \]
M <- matrix(c(2, 0, -1, 0, 0, 0,
5, 0, 0, -1, 0, 0,
0, 1, -3, -2, -2, 0,
0, 1, 0, 0, 0, -1,
0, 3, -4, -4, -1, -2), nrow = 5, byrow = TRUE)
qr(M)$rank
[1] 5
[,1]
[1,] 1
[2,] 40
[3,] 2
[4,] 5
[5,] 12
[6,] 40
Resultado:
\[ (a, b, c, d, e, f) = (1, 40, 2, 5, 12, 40) \]
Equação balanceada:
\[ \text{As}_2\text{S}_5 + 40\,\text{HNO}_3 \rightarrow 2\,\text{H}_3\text{AsO}_4 + 5\,\text{H}_2\text{SO}_4 + 12\,\text{H}_2\text{O} + 40\,\text{NO}_2 \]
A fotossíntese é o processo pelo qual organismos autotróficos convertem dióxido de carbono e água em glicose, liberando oxigênio molecular. A reação completa, com formação e liberação de água, pode ser representada genericamente como:
\[ a\,\text{CO}_2 + b\,\text{H}_2\text{O} \rightarrow c\,\text{C}_6\text{H}_{12}\text{O}_6 + d\,\text{O}_2 + e\,\text{H}_2\text{O} \]
Com as variáveis:
Sistema de equações (elementos C, H, O):
Há 6 incógnitas e 3 equações homogêneas.
# Variáveis: a, b, c, d, e
# Equações:
# 1) a - 6c = 0 (carbono)
# 2) 2b - 12c - 2e = 0 (hidrogênio)
# 3) 2a + b - 6c - 2d - e = 0 (oxigênio)
M <- matrix(c(
1, 0, -6, 0, 0,
0, 2, -12, 0, -2,
2, 1, -6, -2, -1
), nrow = 3, byrow = TRUE)
qr(M)$rank
[1] 3
[,1] [,2]
[1,] -3 0
[2,] -2 -3
[3,] 0 0
[4,] -3 0
[5,] 1 -3
Esse sistema homogêneo não tem solução.
Note que na equação da fotossíntese não balanceada simplificada onde a água aparece apenas como reagente sem contabilizar o excedente de água produzido:
\[ a\,\text{CO}_2 + b\,\text{H}_2\text{O} \rightarrow c\,\text{C}_6\text{H}_{12}\text{O}_6 + d\,\text{O}_2 \]
# Variáveis: a, b, c, d
# Equações:
# 1) a - 6c = 0 (carbono)
# 2) 2b - 12c = 0 (hidrogênio)
# 3) 2a + b - 6c - 2d= 0 (oxigênio)
M <- matrix(c(
1, 0, -6, 0,
0, 2, -12, 0,
2, 1, -6, -2
), nrow = 3, byrow = TRUE)
qr(M)$rank
[1] 3
[,1]
[1,] 6
[2,] 6
[3,] 1
[4,] 6
Sabendo-se que o produto da fotossíntese produz 1 mol de glicose e 6 moles de água. Portanto, \(c=1\) e \(e=6\).
solve a = 6c, 2b = 12c + 2e, 2a + b = 6c + 2d + e, c=1, e=6
Equação balanceada final:
\[ 6\,\text{CO}_2 + 12\,\text{H}_2\text{O} \rightarrow \text{C}_6\text{H}_{12}\text{O}_6 + 6\,\text{O}_2 + 6\,\text{H}_2\text{O} \]
Frequentemente ocorre que os objetos estão relacionados a certas características e que temos que estudar o que chamamos na Seção 3.3 de uma relação. Uma ferramenta útil no estudo das relações é a teoria dos grafos.
Um exemplo bem conhecido são os gráficos químicos ou fórmulas estruturais. Os vértices são átomos e as arestas ligações químicas.
Formamos uma matriz da seguinte maneira: Os átomos são numerados em uma ordem arbitrária. Se houver uma ligação entre dois átomos, atribuímos o número um, caso contrário, o número zero. Os números são coletados em uma matriz quadrada conforme mostrado na Fig. 14.1 (adaptado de Rouvray, 1973). Essas matrizes são chamadas de matrizes de adjacência. A teoria dos grafos é hoje usada em várias áreas, como ecologia, taxonomia, fisiologia e epidemiologia.
Conforme Yoshida (2021, cap. 9), a análise de redes sociais (SNA) pode ser aplicada para entender a conectividade social entre trabalhadores. A conectividade social afeta o estresse de várias maneiras, como hipotetizado por Holt-Lunstad et al. (2010). Embora a conectividade social geralmente seja considerada benéfica para a saúde e bem-estar, nosso estudo inicial mostra que ela também pode ter efeitos negativos, à medida que associamos a conectividade social ao aumento de comportamentos destrutivos, como uso ilegal de drogas e comportamento abusivo. A análise de redes pode ser aplicada à neurologia, genética e redes terroristas.
Fornito et al., 2015, p. 160
Goh et al., 2007, p. 8687
Goh et al., 2007, p. 8687
Em dado país há cinco cidades denominadas A, B, C, D e E que estão interligadas por meio do transporte aéreo de determinado produto. O diagrama de conexão entre as cidades é o seguinte:
arestas <- matrix(c(
"A", "B",
"A", "D",
"B", "C",
"C", "A",
"C", "E",
"D", "C",
"E", "D"
), byrow = TRUE, ncol = 2)
g <- igraph::graph_from_edgelist(arestas, directed = TRUE)
plot(g,
vertex.label.color = "black",
vertex.color = "white",
edge.arrow.size = 0.5,
layout = igraph::layout_with_fr(g))
A B D C E
A 0 1 1 0 0
B 0 0 0 1 0
D 0 0 0 1 0
C 1 0 0 0 1
E 0 0 1 0 0
Qual é o número mínimo de conexões (arestas direcionadas) necessárias para garantir que o grafo seja fortemente conexo, ou seja, que de qualquer cidade seja possível alcançar qualquer outra cidade?
Ou, mais informalmente:
Qual é o menor número de conexões necessárias para garantir que seja possível viajar de qualquer cidade para qualquer outra, respeitando a direção dos voos?
Solução:
A matriz quadrada \(F\) com cinco linhas e cinco colunas contendo os valores 0 e 1 a seguir representa o diagrama de conexões entre as cinco cidades. Por exemplo, o valor 1 na primeira linha e segunda coluna da matriz representa que há a possibilidade de transporte aéreo com origem na cidade A e destino para a cidade B: \(A\to B\). A matriz \(F\) é denominada matriz de adjacência (ou matriz de conectividade direta).
Interpretação da matriz \(F^2\):
O quadrado da matriz de adjacência , isto é, \(F^2\), tem uma interpretação prática importante. Ela informa de quantas maneiras distintas pode-se chegar em cada cidade partindo-se de cada uma delas por meio de exatamente duas conexões. Por exemplo, o valor 2 na primeira linha e terceira coluna da matriz a seguir significa que há exatamente duas maneiras distintas de chegar na cidade C partindo-se da cidade A por meio de duas conexões: \(A\to D\to C\) e \(A\to B\to C\). O valor 1 na segunda linha e primeira coluna significa que há exatamente apenas uma maneira de chegar na cidade A partindo-se de B por meio de duas conexões: \(B\to C\to A\).
Para verificar a acessibilidade entre cidades por múltiplas conexões, computa-se até todos os elementos de \(T^k\) serem não nulos, i.e.:
\[ T^2 = F + F^2 T^3 = F + F^2 + F^3 T^4 = F + F^2 + F^3 + F^4 \]
\[ T^2 =F + F^2= \begin{bmatrix} 0 & 1 & 2 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 2 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ \end{bmatrix} \]
Interpretação da matriz \(T^2\):
A soma da matriz de adjacência com seu quadrado, isto é, \(T^2 = F + F^2\), informa a quantidade de maneiras distintas que se pode chegar a cada cidade partindo-se de cada uma delas por meio de uma ou duas conexões. Por exemplo, o valor 1 na primeira linha e segunda coluna significa que há uma maneira distinta de chegar na cidade B partindo-se da cidade A por meio de uma conexão: \(A\to B\). O valor 1 na segunda linha e primeira coluna significa que há apenas uma maneira de chegar na cidade A partindo-se de B por meio de duas conexões: \(B\to C\to A\).
\[ T^3 = F + F^2 + F^3= \begin{bmatrix} 2 & 1 & 2 & 1 & 2 \\ 1 & 1 & 1 & 2 & 1 \\ 1 & 1 & 3 & 2 & 1 \\ 1 & 1 & 1 & 2 & 1 \\ 1 & 0 & 1 & 1 & 1 \\ \end{bmatrix} \]
\[ T^4 =F + F^2 + F^3 + F^4= \begin{bmatrix} 2 & 3 & 2 & 5 & 2 \\ 1 & 1 & 4 & 2 & 1 \\ 4 & 1 & 3 & 2 & 4 \\ 1 & 1 & 4 & 2 & 1 \\ 1 & 2 & 1 & 3 & 1 \\ \end{bmatrix} \]
sem_zero <- function(M) all(M != 0)
arestas <- matrix(c(
"A", "B",
"A", "D",
"B", "C",
"C", "A",
"C", "E",
"D", "C",
"E", "D"
), byrow = TRUE, ncol = 2)
g <- igraph::graph_from_edgelist(arestas, directed = TRUE)
F <- as.matrix(igraph::as_adjacency_matrix(g,
type = "both",
sparse = FALSE))
F
A B D C E
A 0 1 1 0 0
B 0 0 0 1 0
D 0 0 0 1 0
C 1 0 0 0 1
E 0 0 1 0 0
A B D C E
A 0 1 1 2 0
B 1 0 0 1 1
D 1 0 0 1 1
C 1 1 2 0 1
E 0 0 1 1 0
[1] FALSE
A B D C E
A 2 1 1 2 2
B 1 1 2 1 1
D 1 1 2 1 1
C 1 1 2 3 1
E 1 0 1 1 1
[1] FALSE
A B D C E
A 2 3 5 2 2
B 1 1 2 4 1
D 1 1 2 4 1
C 4 1 2 3 4
E 1 1 3 1 1
[1] TRUE
MatrixPower[{{0,1,0,1,0},{0,0,1,0,0},{1,0,0,0,1}, {0,0,1,0,0},{0,0,0,1,0}},2]
MatrixPower[{{0,1,0,1,0},{0,0,1,0,0},{1,0,0,0,1}, {0,0,1,0,0},{0,0,0,1,0}},3]
MatrixPower[{{0,1,0,1,0},{0,0,1,0,0},{1,0,0,0,1}, {0,0,1,0,0},{0,0,0,1,0}},4]
Essa matriz indica o número de caminhos de comprimento até 4 entre pares de cidades. Como todos os elementos (inclusive fora da diagonal) são positivos, significa que é possível chegar a qualquer cidade a partir de qualquer outra em no máximo 4 conexões.
A diagonal principal indica caminhos que retornam à cidade de origem.
Conclusão: o grafo é conexo em 4 passos.
arestas <- matrix(c(
"A", "B",
"A", "D",
"B", "C",
"C", "A",
"C", "E",
"D", "C",
"E", "D"
), byrow = TRUE, ncol = 2)
g <- igraph::graph_from_edgelist(arestas, directed = TRUE)
plot(g,
vertex.label.color = "black",
vertex.color = "white",
edge.arrow.size = 0.5,
layout = igraph::layout_with_fr(g))
A B D C E
A 0 1 1 0 0
B 0 0 0 1 0
D 0 0 0 1 0
C 1 0 0 0 1
E 0 0 1 0 0
# Função para determinar o menor k tal que T = F + F^2 + ... + F^k não tenha zeros
caminhos_ate_k <- function(F, max_k = 10) {
n <- nrow(F)
T <- F # começa com F^1
P <- F # P armazena a potência atual F^k
for (k in 2:max_k) {
P <- P %*% F # calcula F^k
T <- T + P # soma cumulativa
cat("Passo", k, "- matriz T:\n")
print(T)
if (all(T > 0)) {
return(list(k = k, T = T))
}
}
return(list(k = NA, T = T))
}
# Executar função
resultado <- caminhos_ate_k(F, max_k = 10)
Passo 2 - matriz T:
A B D C E
A 0 1 1 2 0
B 1 0 0 1 1
D 1 0 0 1 1
C 1 1 2 0 1
E 0 0 1 1 0
Passo 3 - matriz T:
A B D C E
A 2 1 1 2 2
B 1 1 2 1 1
D 1 1 2 1 1
C 1 1 2 3 1
E 1 0 1 1 1
Passo 4 - matriz T:
A B D C E
A 2 3 5 2 2
B 1 1 2 4 1
D 1 1 2 4 1
C 4 1 2 3 4
E 1 1 3 1 1
# Imprimir resultado final
if (is.na(resultado$k)) {
cat("Não há k ≤ 10 tal que todas as entradas de T sejam positivas.\n")
} else {
cat("Número de passos para o grafo ser totalmente conexo:", resultado$k, "\n")
}
Número de passos para o grafo ser totalmente conexo: 4
O modelo de insumo-produto pode ser aplicado em problemas de processo de produção industrial. No processo de produção há a fabricação de produtos a partir da matéria-prima e a montagem destes a partir de outros até atingir a produção daqueles demandados pelo mercado. O modelo de insumo-produto aplica-se ao MRP (planejamento das necessidades de material).
Um fabricante vende cinco produtos denotados por A, B, C, D e E.
Os produtos B e E são produtos básicos, fabricados diretamente a partir da matéria-prima, e cada unidade de B ou E consome uma unidade de matéria-prima.
O produto C requer:
O produto A requer:
O produto D requer:
A demanda prevista é:
Os custos unitários de produção dos produtos A, B, C, D e E são, respectivamente:
O custo unitário da matéria-prima utilizada na produção de B e E é igual a 1 unidade monetária, sendo a mesma matéria-prima para ambos.
Passo 1: Construir o diagrama de montagem
O diagrama de montagem a seguir consiste na representação gráfica do esquema de montagem dos produtos.
O círculo representa um produto básico, fabricado diretamente da matéria-prima.
O retângulo representa um produto composto por outro(s) produto(s) composto(s) e/ou básico(s).
O arco orientado (flecha) indica a sequência de montagem e a quantidade de produto empregada.
O nível de montagem representa a complexidade do processo produtivo.
Dessa forma, há três níveis de complexidade no processo produtivo desse problema.
# Nomes dos componentes
componentes <- c("A", "B", "C", "D", "E")
# Matriz de montagem M (linhas: input, colunas: output)
M <- matrix(c(
0, 0, 0, 6, 0, # A
4, 0, 2, 0, 0, # B
1, 0, 0, 2, 0, # C
0, 0, 0, 0, 0, # D
0, 0, 3, 1, 0 # E
), byrow = TRUE, nrow = 5)
M
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 6 0
[2,] 4 0 2 0 0
[3,] 1 0 0 2 0
[4,] 0 0 0 0 0
[5,] 0 0 3 1 0
colnames(M) <- rownames(M) <- componentes
# Gerar arestas ponderadas a partir de M
arestas <- which(M > 0, arr.ind = TRUE)
edges <- data.frame(
from = rownames(M)[arestas[, 1]],
to = colnames(M)[arestas[, 2]],
weight = M[arestas]
)
# Criar grafo direcionado com pesos
g <- igraph::graph_from_data_frame(edges,
directed = TRUE,
vertices = componentes)
# Plotar grafo
plot(
g,
edge.label = igraph::E(g)$weight,
edge.arrow.size = 0.5,
vertex.color = "white",
vertex.label.color = "black",
layout = igraph::layout_as_tree(g,
root = "B",
flip.y = FALSE,
circular = TRUE)
)
A matriz quadrada de montagem \(M\) expressa a relação entre os estágios 0 e 1 do processo produtivo subjacentes ao diagrama de montagem.
O valor \(m_{ij}\) é a quantidade de unidades do produto da linha \(i\) necessária para produzir uma unidade do produto da coluna \(j\), sendo \(i = 1, 2, 3, 4, 5\) e \(j = 1, 2, 3, 4, 5\).
A matriz de estágio 0 é a matriz identidade \(I = M^0\):
\[ I = M^0 = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \]
O estágio produtivo 0 é representado pela matriz identidade \(I\), e cada valor 1 da diagonal principal significa que determinado produto deve ser produzido para satisfazer à sua respectiva demanda por uma unidade dele no mercado. Isto é, a matriz identidade \(I\) expressa a relação entre os estágios 0 e 0 do processo produtivo.
A matriz de montagem total \(T\) representa a quantidade acumulada de produtos (entradas) necessários para a montagem de uma unidade de produto (saída), considerando os estágios indiretos até a profundidade \(M^3\), pois a partir \(M^4=0\):
\[ T = I + M + M^2 + M^3 = \begin{bmatrix} 1 & 0 & 0 & 6 & 0 \\ 6 & 1 & 2 & 40 & 0 \\ 1 & 0 & 1 & 8 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 3 & 0 & 3 & 25 & 1 \\ \end{bmatrix} \]
A interpretação do elemento \(t_{i,j}\) é: quantidade total de unidades do produto da linha \(i\) necessárias para produzir uma unidade do produto da coluna \(j\), considerando toda a cadeia de montagem.
# Definir matriz de montagem M
M <- matrix(c(
0, 0, 0, 6, 0,
4, 0, 2, 0, 0,
1, 0, 0, 2, 0,
0, 0, 0, 0, 0,
0, 0, 3, 1, 0
), nrow = 5, byrow = TRUE)
rownames(M) <- colnames(M) <- c("A", "B", "C", "D", "E")
# Função para calcular potências de M até M^k_max
calcular_potencias_montagem <- function(M, k_max = 4) {
potencias <- list()
potencias[[1]] <- diag(nrow(M)) # M^0 = identidade
cat("M^0: produtos finais diretos (demanda)\n")
print(potencias[[1]])
cat("\n")
for (k in 1:k_max) {
potencias[[k + 1]] <- potencias[[k]] %*% M
cat(paste0("M^", k, ": insumos indiretos de ordem ", k, "\n"))
print(round(potencias[[k + 1]], 4))
cat("\n")
}
return(potencias)
}
# Executar cálculo até M^4
potencias <- calcular_potencias_montagem(M, k_max = 4)
M^0: produtos finais diretos (demanda)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 1
M^1: insumos indiretos de ordem 1
A B C D E
[1,] 0 0 0 6 0
[2,] 4 0 2 0 0
[3,] 1 0 0 2 0
[4,] 0 0 0 0 0
[5,] 0 0 3 1 0
M^2: insumos indiretos de ordem 2
A B C D E
[1,] 0 0 0 0 0
[2,] 2 0 0 28 0
[3,] 0 0 0 6 0
[4,] 0 0 0 0 0
[5,] 3 0 0 6 0
M^3: insumos indiretos de ordem 3
A B C D E
[1,] 0 0 0 0 0
[2,] 0 0 0 12 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 18 0
M^4: insumos indiretos de ordem 4
A B C D E
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
MatrixPower[{{0,0,0,6,0},{4,0,2,0,0},{1,0,0,2,0}, {0,0,0,0,0},{0,0,3,1,0}},2]
MatrixPower[{{0,0,0,6,0},{4,0,2,0,0},{1,0,0,2,0}, {0,0,0,0,0},{0,0,3,1,0}},3]
MatrixPower[{{0,0,0,6,0},{4,0,2,0,0},{1,0,0,2,0}, {0,0,0,0,0},{0,0,3,1,0}},4]
Cálculo da necessidade de matéria-prima:
Uma unidade do produto A necessita de 6 unidades de B e 3 de E.
Uma unidade de C requer 2 unidades de B e 3 de E.
Uma unidade de D consome 40 unidades de B e 25 de E.
Por fim, uma unidade de E necessita de uma unidade de E.
Com base na demanda prevista:
Somando as quantidades de insumos necessários:
Logo, a quantidade total de matéria-prima necessária é 302 + 188 = 490.
A quantidade total de produção necessária para satisfazer à demanda prevista é dada por:
\[ N = T \cdot D = \begin{bmatrix} 46 \\ 302 \\ 59 \\ 5 \\ 188 \\ \end{bmatrix} \]
onde:
\[ D = \begin{bmatrix} 16 \\ 0 \\ 3 \\ 5 \\ 6 \\ \end{bmatrix} \quad \text{correspondente aos produtos} \quad \begin{bmatrix} A \\ B \\ C \\ D \\ E \\ \end{bmatrix} \]
## Determinação de k tal que M^k = 0 e cálculo de N = T × D
# Determinar o menor k tal que M^k = 0
k <- 1
Mk <- M
while (!all(Mk == 0)) {
k <- k + 1
Mk <- expm::'%^%'(M, k)
}
k_min <- k
k_min
[1] 4
# Calcular T = I + M + M^2 + ... + M^{k-1}
T <- diag(5)
for (i in 1:(k_min - 1)) {
T <- T + expm::'%^%'(M, i)
}
colnames(T) <- rownames(T) <- c("A", "B", "C", "D", "E")
T
A B C D E
A 1 0 0 6 0
B 6 1 2 40 0
C 1 0 1 8 0
D 0 0 0 1 0
E 3 0 3 25 1
# Vetor de demanda final D
D <- c(16, 0, 3, 5, 6)
names(D) <- c("A", "B", "C", "D", "E")
# Calcular N = T × D
N <- T %*% D
rownames(N) <- c("A", "B", "C", "D", "E")
N
[,1]
A 46
B 302
C 59
D 5
E 188
Seja:
\[ C = \begin{bmatrix} C_A \\ C_B \\ C_C \\ C_D \\ C_E \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{bmatrix} \]
O custo total de montagem é dado por:
\[ C_M = C^{\prime} N \]
Seja o vetor-coluna dos custos unitários da matéria-prima:
\[ K = \begin{bmatrix} 0 \\ K_B \\ 0 \\ 0 \\ K_E \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 1 \end{bmatrix} \]
O custo total da matéria-prima é:
\[ C_P = K^{\prime} N \]
Portanto, o custo total da produção é:
\[ C_T = C_M + C_P = C^{\prime} N + K^{\prime} N = (C + K)^{\prime} N \]
Substituindo os valores:
\[ \begin{align} C_T &= \left[ C_A,\ C_B + K_B,\ C_C,\ C_D,\ C_E + K_E \right] \begin{bmatrix} 46 \\ 302 \\ 59 \\ 5 \\ 188 \end{bmatrix}\\ &= 46 C_A + 302(C_B + K_B) + 59 C_C + 5 C_D + 188(C_E + K_E)\\ &= 46 \cdot 1 + 302 \cdot 3 + 59 \cdot 3 + 5 \cdot 4 + 188 \cdot 6\\ C_T &= 2277 \end{align} \]
O custo total de produção para satisfazer à demanda é de 2.277 unidades monetárias.
## Cálculo do custo total de produção (incluindo matéria-prima)
# Custos unitários dos produtos A, B, C, D, E
custos_unitarios <- c(A = 1, B = 2, C = 3, D = 4, E = 5)
# Matéria-prima é utilizada na produção de B e E
# Custo da matéria-prima: 1 por unidade
custo_mp_unitario <- 1
quantidade_mp <- N["B", ] + N["E", ]
custo_mp_total <- quantidade_mp * custo_mp_unitario
# Custo de produção dos produtos compostos
custo_produtos <- sum(custos_unitarios * N)
# Custo total
custo_total <- custo_produtos + custo_mp_total
# Resultados
as.numeric(quantidade_mp)
[1] 490
[1] 490
[1] 1787
[1] 2277
\[ T = (I - M)^{-1} \]
A matriz \((I - M)^{-1}\) é chamada matriz de Leontief.
# Matriz M (mesma dos problemas anteriores)
M <- matrix(c(
0, 0, 0, 6, 0,
4, 0, 2, 0, 0,
1, 0, 0, 2, 0,
0, 0, 0, 0, 0,
0, 0, 3, 1, 0
), nrow = 5, byrow = TRUE)
I <- diag(5)
# Matriz de Leontief
T_leontief <- solve(I - M)
T_leontief
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 6 0
[2,] 6 1 2 40 0
[3,] 1 0 1 8 0
[4,] 0 0 0 1 0
[5,] 3 0 3 25 1
# Comparar com T = I + M + M^2 + M^3 + M^4
T_serie <- I
for (k in 1:4) {
T_serie <- T_serie + expm::'%^%'(M, k)
}
# Verificar se são próximas
all.equal(round(T_leontief, 6), round(T_serie, 6))
[1] TRUE
componentes <- c("A", "B", "C", "D", "E")
M <- matrix(c(
0, 0, 0, 6, 0, # A
4, 0, 2, 0, 0, # B
1, 0, 0, 2, 0, # C
0, 0, 0, 0, 0, # D
0, 0, 3, 1, 0 # E
), byrow = TRUE, nrow = 5)
rownames(M) <- colnames(M) <- produtos <- componentes
arestas <- which(M > 0, arr.ind = TRUE)
edge_list <- cbind(produtos[arestas[, 1]], produtos[arestas[, 2]])
pesos <- M[arestas]
g <- igraph::graph_from_edgelist(edge_list, directed = TRUE)
igraph::E(g)$label <- pesos
igraph::E(g)$weight <- pesos
plot(g,
vertex.label.color = "black",
vertex.color = "white",
edge.label = igraph::E(g)$label,
edge.arrow.size = 0.5,
layout = igraph::layout_as_tree(g, root = "B", circular = TRUE),
main = "Diagrama de montagem dos produtos")
A B C D E
A 1 0 0 6 0
B 6 1 2 40 0
C 1 0 1 8 0
D 0 0 0 1 0
E 3 0 3 25 1
[,1]
A 46
B 302
C 59
D 5
E 188
A B C D E
1 2 3 4 5
[1] 2277
Suponha que temos um grafo não direcionado \(G = (V, E)\) com um conjunto de vértices \(V = \{1, 2, \ldots, n\}\). Então a matriz de adjacência de um grafo não direcionado \(G\) é uma matriz \(n \times n\) \(A\) tal que
\[ A = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \]
em que
\[ a_{ij} = \begin{cases} 1 & \text{se há uma aresta não direcionada entre o vértice } i \text{ e o vértice } j \\ 0 & \text{caso contrário} \end{cases} \]
A matriz de adjacência de um grafo não direcionado é uma matriz simétrica.
Fig. 14.1. Matrizes de adjacência de grafo. Se dois vértices estiverem conectados entre si, o número um é atribuído, caso contrário, o número zero.
igraph
packageigraphdata
package# Definir as matrizes de adjacência
adj_matrix1 <- matrix(c(
0, 1, 1, 1, 1,
1, 0, 0, 0, 0,
1, 0, 0, 0, 0,
1, 0, 0, 0, 0,
1, 0, 0, 0, 0
), nrow = 5, byrow = TRUE)
adj_matrix2 <- matrix(c(
0, 1, 0, 0, 0,
1, 0, 1, 0, 0,
0, 1, 0, 1, 0,
0, 0, 1, 0, 1,
0, 0, 0, 1, 0
), nrow = 5, byrow = TRUE)
adj_matrix3 <- matrix(c(
0, 1, 0, 0, 1,
1, 0, 1, 0, 0,
0, 1, 0, 1, 0,
0, 0, 1, 0, 1,
1, 0, 0, 1, 0
), nrow = 5, byrow = TRUE)
g1 <- igraph::graph_from_adjacency_matrix(adj_matrix1, mode = "undirected")
g2 <- igraph::graph_from_adjacency_matrix(adj_matrix2, mode = "undirected")
g3 <- igraph::graph_from_adjacency_matrix(adj_matrix3, mode = "undirected")
plot(g1, main = "Grafo 1", vertex.label = V(g1)$name)
Suponha que temos uma rede \(N = (V, E)\) com um conjunto de vértices \(V = \{1, 2, \ldots, n\}\). Então, a matriz de adjacência de uma rede \(N\) é uma matriz \(n \times n\) \(A\) tal que
\[ A = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \]
onde
\[ a_{ij} = \begin{cases} 1 & \text{se há uma aresta direcionada do vértice } i \text{ para o vértice } j \\ 0 & \text{caso contrário} \end{cases} \]
Para o conjunto de vértices \(V = \{1, 2, 3, 4, 5\}\) e as arestas direcionadas \((1, 2)\), \((1, 3)\), \((1, 4)\), \((4, 3)\), \((3, 2)\), \((2, 5)\), a matriz de adjacência é:
\[ A = \begin{pmatrix} 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \]
A matriz de adjacência não é simétrica.
# Example 160: B
g <- igraph::make_graph(c(1, 2, 1, 3, 1, 4, 4, 3, 3, 2, 2, 5),
directed=TRUE)
plot(g, main = "Grafo Direcionado: Rede")
O grau de um vértice \(v\) em \(V\) para um grafo ou uma rede é o número de arestas adjacentes ao vértice \(v\). A matriz de graus de um grafo ou uma rede com um conjunto de vértices \(V = \{1, 2, \ldots, n\}\) é uma matriz diagonal \(n \times n\) \(D\) tal que
\[ D = \begin{pmatrix} d_1 & 0 & \cdots & 0 & 0 \\ 0 & d_2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & d_{n-1} & 0 \\ 0 & 0 & \cdots & 0 & d_n \end{pmatrix} \]
sendo que \(d_i\) é o grau de um vértice \(i\) em \(V\).
Exemplo A: Grafo triangular não-direcionado
O triângulo com vértices \(V = \{1, 2, 3\}\) mostrado a seguir é um grafo não direcionado.
# Example 158: A
g <- igraph::make_graph(c(1, 2, 2, 3, 3, 1),
directed=FALSE)
plot(g, main = "Grafo Não-Direcionado")
A matriz de graus do triângulo é
\[ D = \begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{pmatrix} \]
[1] 2 2 2
[,1] [,2] [,3]
[1,] 2 0 0
[2,] 0 2 0
[3,] 0 0 2
Exemplo B Grafo direcionado
# Example 160
g <- igraph::make_graph(c(1, 2, 2, 3, 3, 1),
directed=FALSE)
plot(g, main = "Grafo Não-Direcionado")
\[ D = \begin{pmatrix} 3 & 0 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \]
[1] 2 2 2
[,1] [,2] [,3]
[1,] 2 0 0
[2,] 0 2 0
[3,] 0 0 2
A matriz Laplaciana de um grafo ou de uma rede com o conjunto de vértices \(V = \{1, 2, \ldots, n\}\) é uma matriz \(n \times n\) \(L\) tal que \[ L = D - A, \] onde \(D\) é a matriz de graus do grafo ou da rede e \(A\) é a matriz de adjacência do grafo ou da rede.
Exemplo A Do Exemplo 158, a matriz de adjacência do triângulo é
\[ A = \begin{pmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{pmatrix} \]
e a matriz de graus do triângulo é
\[ D = \begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{pmatrix} \]
Portanto, a matriz Laplaciana do triângulo é
\[ L = \begin{pmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{pmatrix} \]
Exemplo B Do Exemplo 160, a matriz de adjacência da rede é
\[ A = \begin{pmatrix} 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \]
e a matriz de graus da rede é
\[ D = \begin{pmatrix} 3 & 0 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \]
Portanto, a matriz Laplaciana da rede é
\[ L = \begin{pmatrix} 3 & -1 & -1 & -1 & 0 \\ 0 & 3 & 0 & 0 & -1 \\ 0 & -1 & 3 & 0 & 0 \\ 0 & 0 & -1 & 2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \]
# Exemplo A: Triângulo
g_triangle <- igraph::make_graph(c(1, 2, 2, 3, 3, 1),
directed=FALSE)
plot(g_triangle, main = "Grafo Não-Direcionado")
A_triangle <- as.matrix(igraph::get.adjacency(g_triangle))
D_triangle <- diag(igraph::degree(g_triangle))
L_triangle <- D_triangle - A_triangle
print("Matriz Laplaciana do Triângulo:")
[1] "Matriz Laplaciana do Triângulo:"
[,1] [,2] [,3]
[1,] 2 -1 -1
[2,] -1 2 -1
[3,] -1 -1 2
[,1] [,2] [,3]
[1,] 2 -1 -1
[2,] -1 2 -1
[3,] -1 -1 2
# Exemplo B: Rede
g_network <- igraph::make_graph(c(1, 2, 1, 3, 1, 4, 4, 3, 3, 2, 2, 5),
directed=TRUE)
plot(g_network, main = "Grafo Direcionado: Rede")
A_network <- as.matrix(igraph::get.adjacency(g_network))
D_network <- diag(igraph::degree(g_network, mode = "out"))
L_network <- D_network - A_network
print("Matriz Laplaciana da Rede:")
[1] "Matriz Laplaciana da Rede:"
[,1] [,2] [,3] [,4] [,5]
[1,] 3 -1 -1 -1 0
[2,] 0 1 0 0 -1
[3,] 0 -1 1 0 0
[4,] 0 0 -1 1 0
[5,] 0 0 0 0 0
[,1] [,2] [,3] [,4] [,5]
[1,] 3 -1 -1 -1 0
[2,] 0 1 0 0 -1
[3,] 0 -1 1 0 0
[4,] 0 0 -1 1 0
[5,] 0 0 0 0 0
A matriz de transição para um grafo não direcionado \(G = (V, E)\) com o conjunto de vértices \(V = \{1, 2, \ldots, n\}\) é uma matriz \(n \times n\) \(P\) tal que
\[ P = \begin{pmatrix} p_{11} & p_{12} & \cdots & p_{1n} \\ p_{21} & p_{22} & \cdots & p_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ p_{n1} & p_{n2} & \cdots & p_{nn} \end{pmatrix} \]
sendo que
\[ p_{ij} = \dfrac{a_{ij}}{d_i} \]
com \(a_{ij}\) sendo o elemento \((i, j)\) da matriz de adjacência do grafo \(G\) e \(d_i\) sendo o grau do vértice \(i\) no grafo \(G\).
Exemplo A: A matriz de transição do triângulo é
\[ P = \begin{pmatrix} 0 & \dfrac{1}{2} & \dfrac{1}{2} \\ \dfrac{1}{2} & 0 & \dfrac{1}{2} \\ \dfrac{1}{2} & \dfrac{1}{2} & 0 \end{pmatrix} \]
A matriz Laplaciana normalizada para um grafo não direcionado \(G = (V, E)\) com o conjunto de vértices \(V = \{1, 2, \ldots, n\}\) é uma matriz \(n \times n\) \(L_P\) tal que
\[ L_P = I_n - P \]
sendo que \(P\) é a matriz de transição do grafo \(G\).
Exemplo A: A matriz de transição do triângulo é
\[ P = \begin{pmatrix} 0 & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & 0 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & 0 \end{pmatrix} \]
Portanto, a matriz Laplaciana normalizada \(L_P\) é
\[ L_P = \begin{pmatrix} 1 & -\frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & 1 & -\frac{1}{2} \\ -\frac{1}{2} & -\frac{1}{2} & 1 \end{pmatrix} \]
g_triangle <- igraph::make_graph(c(1, 2, 2, 3, 3, 1),
directed=FALSE)
plot(g_triangle, main = "Grafo Não-Direcionado")
A_triangle <- as.matrix(igraph::as_adjacency_matrix(g_triangle))
degrees_triangle <- igraph::degree(g_triangle)
P_triangle <- A_triangle / degrees_triangle
I_triangle <- diag(1, nrow = length(degrees_triangle))
L_P_triangle <- I_triangle - P_triangle
print("Matriz de Transição do Triângulo:")
[1] "Matriz de Transição do Triângulo:"
[,1] [,2] [,3]
[1,] 0.0 0.5 0.5
[2,] 0.5 0.0 0.5
[3,] 0.5 0.5 0.0
[1] "Matriz Laplaciana Normalizada do Triângulo:"
[,1] [,2] [,3]
[1,] 1.0 -0.5 -0.5
[2,] -0.5 1.0 -0.5
[3,] -0.5 -0.5 1.0
[,1] [,2] [,3]
[1,] 1.0 -0.5 -0.5
[2,] -0.5 1.0 -0.5
[3,] -0.5 -0.5 1.0
As matrizes Laplacianas nos informam muito sobre as propriedades dos grafos. Para mostrar essas propriedades, precisamos usar ferramentas da álgebra linear. Nesta seção, discutiremos algumas das propriedades bem conhecidas das matrizes Laplacianas. Aqui assumimos que temos grafos não direcionados. Se tivermos um grafo direcionado, ou seja, uma rede, transformamos uma aresta direcionada em uma aresta não direcionada.
Exemplo B Faremos desta rede um grafo não direcionado. Um conjunto de arestas direcionadas consiste em arestas não direcionadas \((1, 2)\), \((1, 3)\), \((1, 4)\), \((4, 3)\), \((3, 2)\), \((2, 5)\). Então, a matriz de adjacência do grafo é
\[ A = \begin{pmatrix} 0 & 1 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \end{pmatrix} \]
e a matriz de graus do grafo é
\[ D = \begin{pmatrix} 3 & 0 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \]
Portanto, a matriz Laplaciana do grafo é
\[ L = \begin{pmatrix} 3 & -1 & -1 & -1 & 0 \\ -1 & 3 & -1 & 0 & -1 \\ -1 & -1 & 3 & -1 & 0 \\ -1 & 0 & -1 & 2 & 0 \\ 0 & -1 & 0 & 0 & 1 \end{pmatrix} \]
Usaremos a matriz Laplaciana do grafo no Exemplo 172 como exemplo contínuo nesta seção.
Teorema 9.3 A matriz Laplaciana \(L\) de um grafo não direcionado \(G\) é simétrica.
Teorema 9.4 A dimensão do espaço nulo da matriz Laplaciana \(L\) de um grafo não direcionado \(G\) é o número de componentes no grafo \(G\).
Usando o Teorema 9.4, podemos ver quantos componentes existem no grafo dado. Isso pode ser especialmente útil quando lidamos com um grafo grande, o que é comum na ciência de dados.
Exemplo B A matriz Laplaciana \(L\) é
\[ L = \begin{pmatrix} 3 & -1 & -1 & -1 & 0 \\ -1 & 3 & -1 & 0 & -1 \\ -1 & -1 & 3 & -1 & 0 \\ -1 & 0 & -1 & 2 & 0 \\ 0 & -1 & 0 & 0 & 1 \end{pmatrix} \]
g <- igraph::make_graph(c(1, 2, 1, 3, 1, 4, 4, 3, 3, 2, 2, 5),
directed=FALSE)
plot(g, main = "Grafo não-Direcionado")
A <- as.matrix(igraph::as_adjacency_matrix(g))
degrees <- igraph::degree(g)
D <- diag(degrees)
L <- D - A
L
[,1] [,2] [,3] [,4] [,5]
[1,] 3 -1 -1 -1 0
[2,] -1 3 -1 0 -1
[3,] -1 -1 3 -1 0
[4,] -1 0 -1 2 0
[5,] 0 -1 0 0 1
[,1] [,2] [,3] [,4] [,5]
[1,] 3 -1 -1 -1 0
[2,] -1 3 -1 0 -1
[3,] -1 -1 3 -1 0
[4,] -1 0 -1 2 0
[5,] 0 -1 0 0 1
[1] "Matriz de Adjacência do Grafo:"
[,1] [,2] [,3] [,4] [,5]
[1,] 0 1 1 1 0
[2,] 1 0 1 0 1
[3,] 1 1 0 1 0
[4,] 1 0 1 0 0
[5,] 0 1 0 0 0
[1] "Matriz de Graus do Grafo:"
[,1] [,2] [,3] [,4] [,5]
[1,] 3 0 0 0 0
[2,] 0 3 0 0 0
[3,] 0 0 3 0 0
[4,] 0 0 0 2 0
[5,] 0 0 0 0 1
[1] "Matriz Laplaciana do Grafo:"
[,1] [,2] [,3] [,4] [,5]
[1,] 3 -1 -1 -1 0
[2,] -1 3 -1 0 -1
[3,] -1 -1 3 -1 0
[4,] -1 0 -1 2 0
[5,] 0 -1 0 0 1
[,1]
[1,] 0.4472136
[2,] 0.4472136
[3,] 0.4472136
[4,] 0.4472136
[5,] 0.4472136
Como há apenas um vetor que gera o espaço nulo de \(L\), a dimensão do espaço nulo de \(L\) é 1.
Teorema 9.5 ([10]) Suponha que \(\lambda_0, \ldots, \lambda_{n-1}\) são autovalores distintos da matriz Laplaciana normalizada \(L_P\) de um grafo não direcionado \(G = (V, E)\) com o conjunto de vértices \(V = \{1, 2, \ldots, n\}\) tal que
\[ \lambda_0 \leq \lambda_1 \leq \ldots \leq \lambda_{n-1} \]
Então,
\[ \sum_{i=0}^{n-1} \lambda_i \leq n \]
e a igualdade vale se e somente se não houver vértices isolados.
Exemplo B A matriz Laplaciana normalizada \(L_P\) é
\[ L_P = \begin{pmatrix} 1 & -\frac{1}{3} & -\frac{1}{3} & -\frac{1}{3} & 0 \\ -\frac{1}{3} & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} \\ -\frac{1}{3} & -\frac{1}{3} & 1 & -\frac{1}{3} & 0 \\ -\frac{1}{2} & 0 & -\frac{1}{2} & 1 & 0 \\ 0 & -1 & 0 & 0 & 1 \end{pmatrix} \]
A soma dos valores próprios da matriz Laplaciana normalizada deste grafo é 5, o que é igual ao número de vértices. Portanto, não há vértice isolado neste grafo.
g <- igraph::make_graph(c(1, 2, 1, 3, 1, 4, 4, 3, 3, 2, 2, 5),
directed=FALSE)
plot(g, main = "Grafo não-Direcionado")
A <- as.matrix(igraph::as_adjacency_matrix(g))
degrees <- igraph::degree(g)
D <- diag(degrees)
L <- D - A
D_inv_sqrt <- diag(1 / sqrt(degrees))
L_P <- D_inv_sqrt %*% L %*% D_inv_sqrt
eigen_values <- eigen(L_P)$values
print("Matriz Laplaciana Normalizada do Grafo:")
[1] "Matriz Laplaciana Normalizada do Grafo:"
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 -0.3333333 -0.3333333 -0.4082483 0.0000000
[2,] -0.3333333 1.0000000 -0.3333333 0.0000000 -0.5773503
[3,] -0.3333333 -0.3333333 1.0000000 -0.4082483 0.0000000
[4,] -0.4082483 0.0000000 -0.4082483 1.0000000 0.0000000
[5,] 0.0000000 -0.5773503 0.0000000 0.0000000 1.0000000
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 -0.3333333 -0.3333333 -0.4082483 0.0000000
[2,] -0.3333333 1.0000000 -0.3333333 0.0000000 -0.5773503
[3,] -0.3333333 -0.3333333 1.0000000 -0.4082483 0.0000000
[4,] -0.4082483 0.0000000 -0.4082483 1.0000000 0.0000000
[5,] 0.0000000 -0.5773503 0.0000000 0.0000000 1.0000000
[1] "Autovalores da Matriz Laplaciana Normalizada:"
[1] 1.767592e+00 1.333333e+00 1.333333e+00 5.657415e-01 2.220446e-15
[1] "Soma dos Autovalores:"
[1] 5
Nesta aplicação prática, usaremos o conjunto de dados “macaque” do
pacote igraphdata
. Este grafo direcionado (rede) é
projetado para modelar as áreas visuotácteis no cérebro do macaco, a fim
de entender como os neurônios no cérebro disparam uns para os outros. O
modelo consiste em 45 áreas (45 vértices na rede) e 463 conexões
direcionadas (463 arestas direcionadas na rede).
Nesta seção, fazemos as seguintes perguntas:
Podemos responder à primeira pergunta usando o Teorema 9.4, que indica que o número de componentes em um grafo não direcionado pode ser determinado pela dimensão do espaço nulo da matriz Laplaciana. Para a segunda pergunta, verificamos se há vértices com grau zero (sem conexões).
library(igraphdata)
data(macaque)
macaque <- igraph::upgrade_graph(macaque)
num_vertices <- igraph::vcount(macaque)
print(paste("Número de vértices no grafo:", num_vertices))
[1] "Número de vértices no grafo: 45"
LM <- igraph::laplacian_matrix(macaque)
LLM <- matrix(LM, LM@Dim[1], LM@Dim[2])
prc <- pracma::nullspace(LLM)
head(prc)
[,1]
[1,] 0.1490712
[2,] 0.1490712
[3,] 0.1490712
[4,] 0.1490712
[5,] 0.1490712
[6,] 0.1490712
Para responder à primeira pergunta, determinamos que há apenas um vetor na base para o espaço nulo da matriz Laplaciana do grafo. Isso significa que há apenas um componente no grafo.
Agora, para responder à segunda pergunta usando o Teorema 9.5, precisamos calcular a matriz Laplaciana normalizada do grafo não direcionado para o conjunto de dados “macaque”.
Ao aplicar o Teorema 9.5, sabemos que a soma dos valores próprios da matriz Laplaciana normalizada deve ser igual ao número de vértices no grafo, se não houver vértices isolados.
Dado que há 45 vértices no conjunto de dados “macaque”, e a soma dos valores próprios da matriz Laplaciana normalizada é 45, podemos concluir que não há vértices isolados no grafo.
LM2 <- igraph::laplacian_matrix(macaque, normalized = TRUE)
LLM2 <- matrix(LM2, LM2@Dim[1], LM2@Dim[2])
EE <- eigen(LLM2)
sum(EE$values)
[1] 45+0i
IGRAPH clustering optimal, groups: 3, mod: 0.38
+ groups:
$`1`
[1] "V1" "V2" "V3" "V3A" "V4t" "VP" "MT" "MSTd/p"
[9] "MSTl" "PO" "LIP" "PIP" "DP" "7a" "FST" "FEF"
$`2`
[1] "V4" "VOT" "PITd" "PITv" "CITd" "CITv" "AITd" "AITv" "STPp" "STPa"
[11] "TF" "TH" "46"
$`3`
[1] "VIP" "3a" "3b" "1" "2" "5" "Ri" "SII" "7b" "4" "6" "SMA"
+ ... omitted several groups/vertices
[1] "Número de vértices no grafo: 45"
igraphdata
: outras aplicações em Biologia# The undirected and connected network of interactions in the immunoglobulin protein. It is made
# up of 1316 vertices representing amino-acids and an edge is drawn between two amino-acids if the
# shortest distance between their C_alpha atoms is smaller than the threshold value θ = 8 Angstrom.
library(igraphdata)
data(immuno)
immuno <- igraph::upgrade_graph(immuno)
num_vertices <- igraph::vcount(immuno)
print(paste("Número de vértices no grafo:", num_vertices))
[1] "Número de vértices no grafo: 1316"
LM <- igraph::laplacian_matrix(immuno)
LLM <- matrix(LM, LM@Dim[1], LM@Dim[2])
prc <- pracma::nullspace(LLM)
head(prc)
[,1]
[1,] 0.02756589
[2,] 0.02756589
[3,] 0.02756589
[4,] 0.02756589
[5,] 0.02756589
[6,] 0.02756589
LM2 <- igraph::laplacian_matrix(immuno, normalized = TRUE)
LLM2 <- matrix(LM2, LM2@Dim[1], LM2@Dim[2])
EE <- eigen(LLM2)
sum(EE$values)
[1] 1316
[1] "Número de vértices no grafo: 1316"
# oc <- igraph::cluster_optimal(immuno)
# print(oc)
# Records of contacts among patients and various types of health care workers in the geriatric unit of a
# hospital in Lyon, France, in 2010, from 1pm on Monday, December 6 to 2pm on Friday, December
# 10. Each of the 75 people in this study consented to wear RFID sensors on small identification
# badges during this period, which made it possible to record when any two of them were in faceto-face contact with each other (i.e., within 1-1.5 m of each other) during a 20-second interval of
# time.
library(igraphdata)
data(rfid)
rfid <- igraph::upgrade_graph(rfid)
num_vertices <- igraph::vcount(rfid)
print(paste("Número de vértices no grafo:", num_vertices))
[1] "Número de vértices no grafo: 75"
LM <- igraph::laplacian_matrix(rfid)
LLM <- matrix(LM, LM@Dim[1], LM@Dim[2])
prc <- pracma::nullspace(LLM)
head(prc)
[,1]
[1,] 0.1154701
[2,] 0.1154701
[3,] 0.1154701
[4,] 0.1154701
[5,] 0.1154701
[6,] 0.1154701
LM2 <- igraph::laplacian_matrix(rfid, normalized = TRUE)
LLM2 <- matrix(LM2, LM2@Dim[1], LM2@Dim[2])
EE <- eigen(LLM2)
sum(EE$values)
[1] 75
[1] "Número de vértices no grafo: 75"
# oc <- igraph::cluster_optimal(rfid)
# print(oc)
# Comprehensive protein-protein interaction maps promise to reveal many aspects of the complex
# regulatory network underlying cellular function.
# This data set was compiled by von Mering et al. (see reference below), combining various sources.
# Only the interactions that have ‘high’ and ‘medium’ confidence are included here.
library(igraphdata)
data(yeast)
yeast <- igraph::upgrade_graph(yeast)
num_vertices <- igraph::vcount(yeast)
print(paste("Número de vértices no grafo:", num_vertices))
[1] "Número de vértices no grafo: 2617"
LM <- igraph::laplacian_matrix(yeast)
LLM <- matrix(LM, LM@Dim[1], LM@Dim[2])
prc <- pracma::nullspace(LLM)
head(prc)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[2,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[3,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[4,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[5,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[6,] 0.001642015 0.0009431883 -0.001996299 0 0.0004257206 -0.001199942
[,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[2,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[3,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[4,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[5,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[6,] 0.000281095 0 0 0 0 -0.000699106 0 0 -0.001417014
[,16] [,17] [,18] [,19] [,20] [,21]
[1,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[2,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[3,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[4,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[5,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[6,] -0.0003481238 0.0003481238 0.001290406 0.002390277 0.00192218 -0.001159913
[,22] [,23] [,24] [,25] [,26]
[1,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[2,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[3,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[4,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[5,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[6,] -0.0004478702 -0.001522364 -0.0002483938 -0.0007427607 0.005981918
[,27] [,28] [,29] [,30] [,31]
[1,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[2,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[3,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[4,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[5,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[6,] -0.003192109 0.001821379 0.004084424 -0.001584059 0.0008077199
[,32] [,33] [,34] [,35] [,36] [,37]
[1,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[2,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[3,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[4,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[5,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[6,] -0.003608892 0.002624766 -0.0001755676 -0.0005091924 0.001496485 0
[,38] [,39] [,40] [,41] [,42] [,43]
[1,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[2,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[3,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[4,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[5,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[6,] 0.003520789 -0.000555508 0 0.002126351 0.0005703229 -0.0006063974
[,44] [,45] [,46] [,47] [,48]
[1,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[2,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[3,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[4,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[5,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[6,] -0.001683246 -0.004235569 0.00188631 7.401527e-05 0.0003023825
[,49] [,50] [,51] [,52] [,53]
[1,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[2,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[3,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[4,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[5,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[6,] 0.0005778596 0.0008226874 0.0009395617 -0.006771942 -0.0003095416
[,54] [,55] [,56] [,57] [,58]
[1,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[2,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[3,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[4,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[5,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[6,] -0.003246146 -0.001943201 -0.0008863114 0.005589436 0.001014074
[,59] [,60] [,61] [,62] [,63] [,64]
[1,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[2,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[3,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[4,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[5,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[6,] 0.003393394 0.002433138 0.002221034 0.001460053 0.001476141 0.003110748
[,65] [,66] [,67] [,68] [,69] [,70]
[1,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[2,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[3,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[4,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[5,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[6,] 0.003736064 0.001101449 0.004573433 0.001079064 0.002865786 0.0009921968
[,71] [,72] [,73] [,74] [,75] [,76]
[1,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[2,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[3,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[4,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[5,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[6,] -0.0001160845 0.00424865 0.004913038 -0.0003380153 0.001735139 0.00306466
[,77] [,78] [,79] [,80] [,81] [,82]
[1,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[2,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[3,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[4,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[5,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[6,] -0.000883367 0.0002783893 -0.0002783893 7.415793e-05 0 -0.00130099
[,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
[1,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[2,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[3,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[4,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[5,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[6,] 0 -0.004181086 0 0 0.002193247 -0.001284092 0 0.0002999713
[,91] [,92]
[1,] -0.0004010922 0
[2,] -0.0004010922 0
[3,] -0.0004010922 0
[4,] -0.0004010922 0
[5,] -0.0004010922 0
[6,] -0.0004010922 0
LM2 <- igraph::laplacian_matrix(yeast, normalized = TRUE)
LLM2 <- matrix(LM2, LM2@Dim[1], LM2@Dim[2])
EE <- eigen(LLM2)
sum(EE$values)
[1] 2617
[1] "Número de vértices no grafo: 2617"
Barabási & Oltvai, 2004, p. 104
# Criar um grafo não direcionado
g1 <- igraph::make_graph(c(1, 2,
2, 3,
4, 5),
directed = FALSE)
summary(g1)
IGRAPH 1b1ca12 U--- 5 3 --
[1] 1
[1] 1 2 1 1 1
[1] 1
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] . 1 . . .
[2,] 1 . 1 . .
[3,] . 1 . . .
[4,] . . . . 1
[5,] . . . 1 .
# Criar um grafo direcionado
g2 <- igraph::make_graph(c(1, 2,
2, 3,
4, 5),
directed = TRUE)
summary(g2)
IGRAPH 1b46ea8 D--- 5 3 --
[1] 2
[1] 1 2 1 1 1
[1] 1
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] . 1 . . .
[2,] . . 1 . .
[3,] . . . . .
[4,] . . . . 1
[5,] . . . . .
# Criar um grafo mais complexo com atributos
g <- igraph::make_graph(~ Alice-Bob:Claire:Frank,
Claire-Alice:Dennis:Frank:Esther,
George-Dennis:Frank,
Dennis-Esther)
g <- igraph::set_vertex_attr(g, "age", value = c(25, 31, 18, 23, 47, 22, 50))
g <- igraph::set_vertex_attr(g, "sex", value = c("f", "m", "f", "m", "m", "f", "m"))
g <- igraph::set_edge_attr(g, "is_formal", value = c(FALSE, FALSE, TRUE, TRUE, TRUE,
FALSE, TRUE, FALSE, FALSE))
summary(g)
IGRAPH 1b5e4b2 UN-- 7 9 --
+ attr: name (v/c), age (v/n), sex (v/c), is_formal (e/l)
[1] 3
Alice Bob Claire Frank Dennis Esther George
3 1 4 3 3 2 2
Claire
4
Claire
3
7 x 7 sparse Matrix of class "dgCMatrix"
Alice Bob Claire Frank Dennis Esther George
Alice . 1 1 1 . . .
Bob 1 . . . . . .
Claire 1 . . 1 1 1 .
Frank 1 . 1 . . . 1
Dennis . . 1 . . 1 1
Esther . . 1 . 1 . .
George . . . 1 1 . .
[1] 6 6 4 3 4 4 4 2 3
IGRAPH clustering optimal, groups: 3, mod: 0.18
+ groups:
$`1`
[1] "Alice" "Bob"
$`2`
[1] "Claire" "Dennis" "Esther"
$`3`
[1] "Frank" "George"
IGRAPH 1b7aad5 UN-- 7 9 --
+ attr: name (v/c)
+ edges from 1b7aad5 (vertex names):
[1] Alice --Bob Alice --Claire Alice --Frank Claire--Frank Claire--Dennis
[6] Claire--Esther Frank --George Dennis--Esther Dennis--George
Na Seção 11.3, estudamos um processo determinístico de nascimento e morte. Neste modelo não consideramos o fato de que as taxas de natalidade e mortalidade são dependentes da idade dos indivíduos.
Lewis (1942) e Leslie (1945) introduziram um modelo determinístico que leva em consideração a estrutura etária da população.
Para simplificar o tratamento matemático, consideramos o processo de nascimento e morte em etapas de duração constante.
Seja \(\Delta t\) um intervalo de tempo adequadamente escolhido. Para uma população humana podemos escolher \(\Delta t=5\) anos ou, para uma análise mais precisa, \(\Delta t=1\) ano. Com o valor escolhido de \(\Delta t\) consideramos a estrutura etária da população nos momentos \(t = 0, \Delta t, 2\Delta t, \ldots\).
Introduzimos também as faixas etárias \(x = 0, 1, 2, \ldots,m\): grupo \(x = 0\) contém as idades de \(0\) a \(\Delta t\), grupo \(x = 1\) as idades de \(\Delta t\) a \(2\Delta t\), grupo \(x = 2\) as idades de \(2\Delta t\) a \(3\Delta t\) etc. A última faixa etária possível é denotada por \(x = m\).
Observe que o intervalo de grupo tem o mesmo comprimento que o intervalo entre instantes de tempo consecutivos da população.
Em uma população com dois sexos, precisamos considerar apenas as fêmeas. A estrutura etária dos machos é de menor importância.
No momento \(t=k\Delta t\), \(k=0,1,2,\ldots\), o tamanho da população feminina é representado pelo vetor coluna \(n_i=\left[n_{1k}, n_{2k}, \cdots, n_{mk}\right]^{\prime}\). Portanto, \(n_{xk}\) é o número de fêmeas do grupo etário \(x\) no momento \(t=k\Delta t\).
\[ n_0 = \begin{pmatrix} n_{00} \\ n_{10} \\ n_{20} \\ \vdots \\ n_{m0} \end{pmatrix}, \quad n_1 = \begin{pmatrix} n_{01} \\ n_{11} \\ n_{21} \\ \vdots \\ n_{m1} \end{pmatrix}, \quad n_2 = \begin{pmatrix} n_{02} \\ n_{12} \\ n_{22} \\ \vdots \\ n_{m2} \end{pmatrix}, \quad \ldots \]
Seja \(F_0, F_1, F_2 , \ldots, F_m\) o número médio de filhas nascidas no intervalo de tempo \(\Delta t\) para uma mulher de idade \(0, 1,2, \ldots, m\), respectivamente.
As fêmeas podem não ser reprodutivas em todas as idades. Pode haver uma fase pré-reprodutiva e uma fase pós-reprodutiva. Assim, alguns dos \(F_x\) , digamos \(F_0\), \(F_1\) e \(F_m\), podem ser zero. O número total de filhas nascidas durante o primeiro intervalo de tempo de \(t = 0\) a \(t = \Delta t\) é então \(\sum_{x=0}^{m}{F_x n_{x0}}\).
Agora admitimos que por \(F_0, F_1, F_2, ..., F_m\) contamos somente as meninas que sobreviveram até que o intervalo de tempo no qual elas nasceram tivesse passado. Ao final desse intervalo todas serão consideradas como de idade \(x=0\). Assim, o número \(n_{01}\) de fêmeas do grupo \(x=0\) no tempo \(t = \Delta t\) é igual ao resultado da fórmula \(\sum_{x=0}^{m}{F_x n_{x0}}\), i.e.:
\[ \sum_{x=0}^{m}{F_x n_{x0}}=n_{01} \]
Por \(P_x\) denotamos a probabilidade de uma fêmea da faixa etária \(x\) sobreviver e entrar na faixa etária \(x + 1\). Assim, para \(x = 0,1, \ldots, m -1\) obtemos
\[ P_xn_{xk}=n_{x+1,k+1} \]
Observe que \(P_m=0\).
No caso especial \(k=0\) obtemos \(P_xn_{x0}=n_{x+1,1}\).
A transição da população de \(t = 0\) para \(t = \Delta t\) pode ser resumida pela multiplicação de matrizes como segue
\[ \begin{pmatrix} F_0 & F_1 & \cdots & F_{m-1} & F_m \\ P_0 & 0 & \cdots & 0 & 0 \\ 0 & P_1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & P_{m-1} & 0 \end{pmatrix} \begin{pmatrix} n_{00} \\ n_{10} \\ n_{20} \\ \vdots \\ n_{m0} \end{pmatrix} = \begin{pmatrix} n_{01} \\ n_{11} \\ n_{21} \\ \vdots \\ n_{m1} \end{pmatrix} \]
Ou em notação matricial:
\[ \large Mn_{0}=n_{1} \]
Supondo que \(F_x\) e \(P_x\) permaneçam constantes com o passar do tempo, podemos repetir o mesmo procedimento. Portanto, \(Mn_1=n_2, Mn_2=n_3, \ldots\).
Então \(n_2=M(Mn_0)=M^2n_0\). Por indução,
\[ \large n_g=M^gn_0 \]
Dado o tamanho inicial da população e a estrutura etária e dada a chamada matriz de projeção \(M\), podemos calcular o tamanho e a estrutura etária da população futura por \(n_g=M^gn_0\).
Dado \(n_0\) apenas, podemos perguntar:
Que propriedade \(M\) deve ter para que a população seja estável em tamanho, isto é, se \(n_0 = n_1 = n_2 = \cdots\)?
Para que a população seja estável em tamanho é necessário que o vetor populacional \(n_0\) seja autovetor da matriz de projeção \(M\) associado ao autovalor \(\lambda = 1\).
Matematicamente, a condição é:
\[ M n_0 = \lambda n_0 \quad \text{com} \quad \lambda = 1. \]
Se \(\lambda = 1\), a população mantém o mesmo tamanho total em cada geração. Se \(\lambda > 1\), a população cresce exponencialmente. Se \(\lambda < 1\), a população decresce ao longo do tempo.
Portanto, a matriz \(M\) deve ter autovalor dominante igual a 1.
Também podemos estar interessados em uma distribuição etária estável. Por essa noção, queremos dizer proporções constantes entre diferentes faixas etárias. Esta questão leva à equação
\[ n_{g+1}=\lambda n_g \]
ou
\[ Mn_g=\lambda n_g \]
Sendo que \(\lambda\) denota um número positivo que explica o possível crescimento populacional. Se nos for dado \(M\), a solução de \(Mn_g=\lambda n_g\) para \(n_g\) e \(\lambda\) requer ferramentas avançadas de teoria de matrizes, como o conceito de equação característica de uma matriz (consulte a Seção 14.9).
# Delta t = 1 ano
# Matriz de Leslie
M <- matrix(c(
0, 2, 1,
0.7, 0, 0,
0, 0.9, 0
), nrow = 3, byrow = TRUE)
# Vetor populacional inicial n0
n0 <- c(100, 50, 20)
# Número de períodos a simular
g <- 10
# Armazenar resultados
n_matrix <- matrix(0, nrow = 3, ncol = g + 1)
n_matrix[, 1] <- n0
# Calcular n_k = M^k * n0 para k = 1, 2, ..., g
for (k in 1:g) {
Mk <- expm::'%^%'(M, k)
n_matrix[, k + 1] <- Mk %*% n0
}
# Nomear linhas e colunas
rownames(n_matrix) <- paste0("faixa", 0:2)
colnames(n_matrix) <- paste0("t", 0:g)
# Visualizar
print(round(n_matrix))
t0 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10
faixa0 100 120 185 231 335 440 614 827 1137 1544 2112
faixa1 50 70 84 130 162 234 308 430 579 796 1081
faixa2 20 45 63 76 117 146 211 277 387 521 716
O exemplo a seguir é adaptado de Thrall et al. (1967, PL 5).
É bem conhecido que poluentes venenosos (como DDT ou mercúrio) se acumulam nas cadeias alimentares.
Escolhemos estudar um sistema ecológico particular com os seguintes três elos de uma cadeia:
Podemos perguntar:
Qual é a quantidade de planta \(p_i\) que é ingerida indiretamente pelo carnívoro \(c_j\) durante uma determinada estação?
Para responder a esta questão, introduzimos a seguinte matriz \(X\) para a transição da conexão 1 para a 2:
\[ \begin{array}{c|cccc} & a_1 & a_2 & \cdots & a_s \\ \hline p_1 & x_{11} & x_{12} & \cdots & x_{1s} \\ p_2 & x_{21} & x_{22} & \cdots & x_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ p_r & x_{r1} & x_{r2} & \cdots & x_{rs} \end{array} = X \]
Aqui \(x_{11}\) denota a quantidade média (em g ou kg) de planta \(p_1\) consumida por cada indivíduo da espécie \(a_1\) durante a estação.
Geralmente, \(x_{ik}\) é a quantidade média de planta \(p_i\) consumida por cada indivíduo da espécie \(a_k\).
Também definimos uma matriz \(Y\) para a transição da conexão 2 para a 3:
\[ \begin{array}{c|cccc} & c_1 & c_2 & \cdots & c_t \\ \hline a_1 & y_{11} & y_{12} & \cdots & y_{1t} \\ a_2 & y_{21} & y_{22} & \cdots & y_{2t} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_s & y_{s1} & y_{s2} & \cdots & y_{st} \end{array} = Y \]
Aqui \(y_{11}\) denota o número de animais da espécie \(a_1\) devorados por todos os indivíduos da espécie \(c_1\) juntos. Geralmente, \(y_{kj}\) é o número de animais da espécie \(a_k\) devorados por carnívoros da espécie \(c_j\) durante a estação. Observe que \(y_{kj}\) é uma variável discreta (contagem), enquanto \(x_{ik}\) é uma variável contínua (real).
Considere agora os animais da espécie \(c_1\). Alimentando-se de espécies \(a_1\), eles consomem indiretamente a quantidade \(x_{11}y_{11}\) da planta \(p_1\). Alimentando-se da espécie \(a_2\), eles consomem \(x_{12}y_{21}\) da planta \(p_1\) etc. A quantidade total média de planta \(p_1\) consumida indiretamente por todos os carnívoros da espécie \(c_1\) é, portanto, \(\sum_{i=1}^{s}{x_{1i}y_{i1}}\).
O resultado é um produto interno, mais precisamente, a primeira linha de \(X\) multiplicada pela primeira coluna de \(Y\).
O resultado pode ser generalizado rapidamente. A quantidade de planta \(p_i\) consumida indiretamente pelo carnívoro \(c_j\) é o produto da \(i\)-ésima linha de \(X\) e da \(j\)-ésima coluna de \(Y\). Obtemos todos os resultados particulares do produto da matriz
\[ XY \]
Isso responde a nossa pergunta.
Considere:
A matriz \(X\) representa o consumo de plantas por herbívoros.
A matriz \(Y\) representa quantos herbívoros são comidos por carnívoros.
O vetor \(z\) contém a quantidade média de mercúrio por unidade de planta.
Queremos calcular a quantidade total de mercúrio ingerida por carnívoros, planta a planta. Para isso:
sweep()
sobre
as linhas da matriz \(XY\)# Matriz X: consumo de plantas por herbívoros
X <- matrix(c(
2, 1, # p1 por a1, a2
0, 3, # p2 por a1, a2
1, 2 # p3 por a1, a2
), nrow = 3, byrow = TRUE)
rownames(X) <- c("p1", "p2", "p3")
colnames(X) <- c("a1", "a2")
# Matriz Y: herbívoros comidos por carnívoros
Y <- matrix(c(
5, 1, # a1 por c1, c2
3, 2 # a2 por c1, c2
), nrow = 2, byrow = TRUE)
rownames(Y) <- c("a1", "a2")
colnames(Y) <- c("c1", "c2")
# Vetor z: contaminação (mg) por unidade de planta
z <- c(0.2, 0.5, 0.1)
names(z) <- rownames(X)
# Produto XY: total indireto de planta consumida por carnívoros
XY <- X %*% Y
# Aplicar vetor de contaminação com sweep (linha por linha)
# Hg[i,j] = z[i] * XY[i,j]
Hg <- sweep(XY, MARGIN = 1, STATS = z, FUN = "*")
# Nomear matriz
rownames(Hg) <- rownames(X)
colnames(Hg) <- colnames(Y)
# Resultado
print(round(Hg, 2))
c1 c2
p1 2.6 0.8
p2 4.5 3.0
p3 1.1 0.5
Considere um locus bialélico com alelos \(A\) e \(a\), com frequências gênicas \(p = P(A)\) e \(q = 1 - p = P(a)\) em uma população panmítica.
Sejam \(R_1\) e \(R_2\) dois indivíduos relacionados por transmissão direta de um alelo idêntico por descendência (IBD). Isso inclui relações como pai-filho, mãe-filho ou avô-neto. O objetivo é modelar a probabilidade condicional do genótipo de \(R_2\) dado o genótipo de \(R_1\), assumindo que compartilham um único alelo IBD.
Esse sistema pode ser representado como uma cadeia de Markov de primeira ordem, com estados discretos \(\{AA, Aa, aa\}\), onde a matriz de transição \(T\) expressa a distribuição condicional dos genótipos de \(R_2\) dado o genótipo de \(R_1\), sob a restrição de compartilhamento de um alelo IBD:
\[ \begin{array}{c|ccc} \text{Genótipo de } R_1 \backslash R_2 & AA & Aa & aa \\ \hline AA & p & q & 0 \\ Aa & \frac{1}{2}p & \frac{1}{2} & \frac{1}{2}q \\ aa & 0 & p & q \end{array} = T \]
\[ T = \begin{pmatrix} P(AA|AA) & P(Aa|AA) & P(aa|AA) \\ P(AA|Aa) & P(Aa|Aa) & P(aa|Aa) \\ P(AA|aa) & P(Aa|aa) & P(aa|aa) \end{pmatrix} = \begin{pmatrix} p & q & 0 \\ \frac{1}{2}p & \frac{1}{2} & \frac{1}{2}q \\ 0 & p & q \end{pmatrix} \]
Indivíduo ----- Pool genético R1 ----- P(A)=p P(a)=q | R2
P(A1)=1/2 P(A2)=1/2 A1A2 ----- P(A)=p P(a)=q | A1A: 1/2 p + A2A: 1/2 p = p
Aa ----- P(A)=p P(a)=q | Aa: 1/2 p + 1/2 q = 1/2
Resultado: P(AA|AA) = 1 × p = p. Com certeza um alelo A é oriundo do genitor \(R_1\). Com probabilidade p o outro alelo A é oriundo do outro genitor.
Se \(R_1\) for do genótipo AA, então \(R_2\) deve ter o alelo A em comum com \(R_1\). Como o segundo alelo é independente do alelo conhecido, esse segundo alelo é A com probabilidade \(p\) ou a com probabilidade \(q\). O parente \(R_2\) não pode ser do genótipo aa. Esses argumentos provam a primeira linha de \(T\).
Se \(R_1\) é do genótipo Aa, então \(R_2\) deve ter o alelo A ou o alelo a em comum com \(R_1\), ambos com probabilidade \(1/2\). Como o segundo alelo é independente do primeiro alelo, esse segundo alelo é A com probabilidade \(p\) ou a com probabilidade \(q\). Portanto, \(R_2\) é do genótipo AA com probabilidade \(\frac{1}{2}p\), do genótipo Aa com probabilidade \(\frac{1}{2}p+\frac{1}{2}q=\frac{1}{2}\) ou do genótipo aa com probabilidade \(\frac{1}{2}q\). Esses argumentos provam a segunda linha de \(T\). A prova da terceira linha é bastante análoga.
A construção da matriz é baseada na suposição de que o segundo alelo de \(R_2\) é amostrado independentemente da população, enquanto o primeiro é herdado diretamente de \(R_1\) (alelo IBD).
Essa matriz é estocástica por linhas e define uma cadeia de Markov com transições entre genótipos ao longo de gerações. O produto de potências da matriz, \(T^n\), fornece a distribuição dos genótipos em descendentes de \(n\)-ésima geração, assumindo herança mendeliana e transmissão de um único alelo IBD.
Exemplo: para um avô \(R_1\) do genótipo \(AA\), o vetor de probabilidades para os genótipos dos netos \(R_3\) é dado por
\[ \text{Linha}_\text{AA} \cdot T^2 = (p, q, 0) \cdot T^2 \]
Mais geralmente, o sistema converge para uma distribuição estacionária \(\pi\), que satisfaz \(\pi = \pi T\), caracterizando a distribuição limite dos genótipos em descendentes sob transmissão repetida de um alelo IBD e recombinação populacional no segundo alelo.
Com \(p = 0.7\), o vetor estacionário é:
\[ \pi = (0.49,\ 0.42,\ 0.09) \]
ou seja, após várias gerações, os genótipos \(AA\), \(Aa\) e \(aa\) estabilizam nessas proporções, refletindo um equilíbrio entre transmissão IBD e amostragem aleatória do alelo não IBD.
Com \(P(A) = 0.7\) e \(P(a) = 0.3\), a tabela de Punnett combina os alelos parentais:
A (0.7) | a (0.3) | |
---|---|---|
A (0.7) | AA (0.49) | Aa (0.21) |
a (0.3) | Aa (0.21) | aa (0.09) |
Somando os genótipos:
Portanto, a distribuição genotípica \((0.49, 0.42, 0.09)\) da tabela de Punnett coincide com o vetor estacionário e com o equilíbrio de Hardy-Weinberg.
Considere um determinado locus gênico com possíveis alelos A e a. Suponha que em uma população o alelo A esteja presente com probabilidade \(p\) e o alelo a com probabilidade \(q= 1-p\).
{{p,q,0},{p/2,1/2,q/2},{0,p,q}}
{{p,1-p,0},{p/2,1/2,(1-p)/2},{0,p,1-p}}
Agora damos um passo adiante e investigamos a relação entre avós e netos. Suponha, por exemplo, que um avô seja do genótipo AA. Então derivamos facilmente as probabilidades para os diferentes genótipos para os netos pelo mesmo método:
Então, para os netos
\[ \begin{align} P(AA) &= pp + q\dfrac{1}{2}p+00=p^2+\dfrac{1}{2}pq\\ P(Aa) &= pq + q\dfrac{1}{2}+0p=pq + \dfrac{1}{2}p\\ P(aa) &= p0 + q\dfrac{1}{2}q+0q= \dfrac{1}{2}q^2 \end{align} \]
{p,q,0} . {{p,q,0},{p/2,1/2,q/2},{0,p,q}}
Essas expressões podem ser interpretadas como “primeira linha de \(T\) multiplicada pela primeira, segunda e terceira coluna de \(T\), respectivamente”. Em geral, obtemos as probabilidades dos genótipos AA, Aa, aa para netos simplesmente por \(TT= T^2\).
O resultado pode ser escrito na forma:
{{p,q,0},{p/2,1/2,q/2},{0,p,q}}^2
{{p,1-p,0},{p/2,1/2,(1-p)/2},{0,p,1-p}}^2
Seguindo para bisnetos teríamos que elaborar a multiplicação da matriz considerando \(TT^2= T^3\)
{{p,q,0},{p/2,1/2,q/2},{0,p,q}}^3
{{p,1-p,0},{p/2,1/2,(1-p)/2},{0,p,1-p}}^3
{{p,1-p,0},{p/2,1/2,(1-p)/2},{0,p,1-p}}^6
{{p,1-p,0},{p/2,1/2,(1-p)/2},{0,p,1-p}}^6 /. p=0.7
Solução Analítica do Vetor Estacionário
O vetor estacionário \(\pi\) é o autovetor associado ao autovalor 1 da matriz de transição \(T\).
Dada a matriz de transição \(T\):
\[ T = \begin{pmatrix} p & q & 0 \\ \frac{1}{2} p & \frac{1}{2} & \frac{1}{2} q \\ 0 & p & q \end{pmatrix} \]
O vetor estacionário \(\pi = (\pi_1, \pi_2, \pi_3)\) satisfaz a equação de equilíbrio:
\[ \pi T = \pi \]
A equação \(\pi T = \pi\) pode ser reescrita como:
\[ \pi T = \pi \cdot 1 \]
Isso significa que \(\pi\) é um autovetor da matriz \(T\) associado ao autovalor \(\lambda = 1\).
Além disso, a equação \(\pi = \pi T\) implica que para qualquer \(n\):
\[ \pi T^n = \pi \]
Isso significa que a distribuição estacionária \(\pi\) permanece invariante sob múltiplas aplicações da matriz de transição \(T\).
As equações de equilíbrio para o vetor estacionário \(\pi = (\pi_1, \pi_2, \pi_3)\) são:
\[ \pi = \pi T \]
Isso resulta nas seguintes equações de equilíbrio:
Resolvendo este sistema, obtemos:
\[ \pi_1 = \frac{p }{2 q}\pi_2 \]
\[ \pi_3 = \frac{q }{2 p} \pi_2 \]
\[ \pi_2 = \frac{2 pq}{1-q^2} \]
Substituindo os valores \(p = 0.7\) e \(q = 0.3\), obtemos diretamente a distribuição estacionária:
\[ \pi = (0.49, 0.42, 0.09) \]
Essa solução analítica mostra que, após um número suficiente de passos, sistema estará 49% no estado \(AA\), 42% em \(Aa\), e 9% em \(aa\).
Além disso, com o chamado método ITO, é fácil obter resultados correspondentes para relacionamentos como irmão-irmão, irmão-meio-irmão, tio-sobrinho, primo-primo (cr. Problema 14.3.3).
# Parâmetros
p <- 0.7
q <- 1 - p
# Matriz de transição T
T <- matrix(c(
p, q, 0,
p/2, 1/2, q/2,
0, p, q
), nrow = 3, byrow = TRUE)
# Vetor inicial
v <- c(p, q, 0)
# Multiplicação vetor linha por T
resultado1 <- v %*% T
print(resultado1)
[,1] [,2] [,3]
[1,] 0.595 0.36 0.045
[,1] [,2] [,3]
[1,] 0.595 0.36 0.045
[2,] 0.420 0.46 0.120
[3,] 0.245 0.56 0.195
[,1] [,2] [,3]
[1,] 0.5425 0.39 0.0675
[2,] 0.4550 0.44 0.1050
[3,] 0.3675 0.49 0.1425
[,1] [,2] [,3]
[1,] 0.4965625 0.41625 0.0871875
[2,] 0.4856250 0.42250 0.0918750
[3,] 0.4746875 0.42875 0.0965625
[,1] [,2] [,3]
[1,] 0.4965625 0.41625 0.0871875
[2,] 0.4856250 0.42250 0.0918750
[3,] 0.4746875 0.42875 0.0965625
# Vetor estacionário: autovalor 1
eig <- eigen(t(T)) # transposta pois queremos solução à esquerda: π T = π
index <- which(abs(eig$values - 1) < 1e-8)
pi_raw <- Re(eig$vectors[, index])
pi_stationary <- proportions(pi_raw) # normaliza
print(pi_stationary)
[1] 0.49 0.42 0.09
Vamos construir a matriz de transição \(T'\) para o caso irmão–irmão, usando o método ITO. Assumimos herança mendeliana e que os irmãos compartilham em média:
Esses coeficientes de identidade determinam a probabilidade condicional do genótipo de um irmão dado o genótipo do outro.
Seja \(T^{\prime}\) a matriz com linhas representando o genótipo de \(R_1\) (irmão 1) e colunas o de \(R_2\) (irmão 2). Ela será uma média ponderada:
\[ T^{\prime} = \dfrac{1}{4} M_2 + \dfrac{1}{2} M_1 + \dfrac{1}{4} M_0 \]
onde:
As matrizes:
\[ M_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad M_0 = \begin{pmatrix} p^2 & 2pq & q^2 \\ p^2 & 2pq & q^2 \\ p^2 & 2pq & q^2 \end{pmatrix}, \quad M_1 = \begin{pmatrix} p & q & 0 \\ \dfrac{1}{2}p & \dfrac{1}{2} & \dfrac{1}{2}q \\ 0 & p & q \end{pmatrix} = T \text{ do caso pai-filho} \]
Logo, a matriz final:
\[ T^{\prime} = \dfrac{1}{4} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + \dfrac{1}{2} \begin{pmatrix} p & q & 0 \\ \dfrac{1}{2}p & \dfrac{1}{2} & \dfrac{1}{2}q \\ 0 & p & q \end{pmatrix} + \dfrac{1}{4} \begin{pmatrix} p^2 & 2pq & q^2 \\ p^2 & 2pq & q^2 \\ p^2 & 2pq & q^2 \end{pmatrix} \]
A matriz de transição \(T'\) para o caso irmão–irmão com \(p = 0.7\) e \(q = 0.3\) é:
\[ T^{\prime} = \begin{pmatrix} 0.7225 & 0.255 & 0.0225 \\ 0.2975 & 0.605 & 0.0975 \\ 0.1225 & 0.455 & 0.4225 \end{pmatrix} \]
O vetor estacionário correspondente é:
\[ \pi = (0.49,\ 0.42,\ 0.09) \]
Ou seja, mesmo com uma estrutura de transição diferente (irmão–irmão vs. pai–filho), a distribuição estacionária dos genótipos é a mesma. Isso mostra que o equilíbrio de Hardy-Weinberg emerge como distribuição limite sob múltiplas gerações, independentemente do tipo de relação de parentesco entre pares consecutivos, desde que haja aleatoriedade suficiente no alelo não IBD.
Uma das moléculas de RNA ficou famosa por ser capaz de se replicar em um tubo de ensaio sem a ajuda de células vivas (Mills et al., 1973). Como outras moléculas de RNA, é uma longa cadeia de apenas quatro constituintes, as chamadas bases: A = adenina, C = citosina, G = guanina, U = uracila.
Fig. 14.2. Uma molécula de RNA, 218 nucleotídeos de comprimento. As quatro bases são A (adenina), C (citosina), G (guanina), U (uracila). Como mostrado no texto, a transição de nucleotídeos consecutivos não pode ser explicada por um processo de Markov.
A molécula é representada na Fig. 14.2. À primeira vista, reconhecemos regularidades na sequência de bases no início (à esquerda) e no final (à direita).
Caso contrário, as bases parecem seguir umas às outras em uma ordem “aleatória”.
Para investigar o tipo de aleatoriedade podemos usar uma técnica que foi inventada por A. A. Markov. Consideramos a cadeia como resultado de um processo estocástico onde cada elo da cadeia está em um dos quatro estados A, C, G, U. A mudança de um link para o seguinte, ou seja, de um estado para outro, é chamada de transição. Com quatro estados temos 16 transições diferentes \(A \to A\), \(A \to C\), \(A \to G\), \(A \to U\), \(C \to A\), \(C \to C\) etc.
As 16 transições possíveis ocorrem com as seguintes frequências:
\[ \begin{array}{c|cccc|c} & A & C & G & U & \text{total} \\ \hline A \rightarrow & 4 & 14 & 13 & 0 & 31 \\ C \rightarrow & 7 & 23 & 32 & 11 & 73 \\ G \rightarrow & 18 & 25 & 26 & 12 & 81 \\ U \rightarrow & 2 & 12 & 9 & 9 & 32 \\ \end{array} \]
As transições que ocorrem com mais frequência são \(C \to G\), \(G \to G\), \(G \to C\) e \(C \to C\), enquanto \(A \to A\) e \(U \to A\) são raras. A transição \(A \to U\) nunca ocorre em nossa sequência especial de RNA.
Das frequências (absolutas) derivamos as frequências relativas em cada linha:
\[ \begin{array}{c|cccc|c} & A & C & G & U & \text{total} \\ \hline A \rightarrow & 0.129 & 0.452 & 0.419 & 0.000 & 1.000 \\ C \rightarrow & 0.096 & 0.315 & 0.438 & 0.151 & 1.000 \\ G \rightarrow & 0.222 & 0.309 & 0.321 & 0.148 & 1.000 \\ U \rightarrow & 0.063 & 0.375 & 0.281 & 0.281 & 1.000 \\ \end{array} \]
\[ F=\left[\begin{array}{rrrr} 0.129 & 0.452 & 0.419 & 0.000\\ 0.096 & 0.315 & 0.438 & 0.151\\ 0.222 & 0.309 & 0.321 & 0.148\\ 0.063 & 0.375 & 0.281 & 0.281 \end{array}\right] \]
Observe que a soma dos elementos em cada linha é igual a 1.
Podemos teorizar e pensar em um processo químico que gera tais transições com certas probabilidades, as probabilidades de transição.
Eles formam uma matriz de transição teórica:
\[ P=\left[\begin{array}{rrrr} p_{11} & p_{12} & p_{13} & p_{14}\\ p_{21} & p_{22} & p_{23} & p_{24}\\ p_{31} & p_{32} & p_{33} & p_{34}\\ p_{41} & p_{42} & p_{43} & p_{44} \end{array}\right] \]
\(p_{11}\) é a probabilidade de transição \(A \to A\) e é estimada pelo valor empírico 0.129. Da mesma forma, \(p_{12}\) é a probabilidade de transição \(A \to C\) e é estimado pelo valor empírico 0.452 etc. Para cada linha postulamos que:
\[ \sum_{k=1}^{4}{p_{ik}}=1, \; i=1,2,3,4 \]
Agora estendemos nossa investigação para transições de duas etapas, como \(A \to C \to A\), \(G \to C \to U\) etc. e podemos perguntar:
Quais são as frequências relativas e/ou probabilidades das 16 transições de duas etapas \(A \to \cdots \to A\), \(A \to \cdots \to C\) etc?
O ponto significa que não estamos interessados no estado específico intermediário. A transição \(A \to \cdots \to A\) na verdade contém os quatro casos a seguir:
\[ A\to A\to A\\ A\to C\to A\\ A\to G\to A\\ A\to U\to A \]
Em cadeias de Markov, assumimos que as transições consecutivas são independentes umas das outras no sentido da teoria da probabilidade (Seção 13.6). Esta é a propriedade de Markov. Em uma escala de tempo, a propriedade de Markov significa que os eventos futuros dependem apenas do tempo presente, não do passado.
Sob a suposição de cadeias de Markov a regra de multiplicação vale, e obtemos:
\[ P(A\to A\to A) = P(A\to A)\times P(A\to A) = p_{11}p_{11}\\ P(A \to C \to A) = P(A \to C)\times P(C \to A) = p_{12}p_{21}\\ P(A \to G\to A) = P(A \to G)\times P(G\to A) = p_{13}p_{31}\\ P(A \to U \to A) = P(A \to U) \times P(U\to A) = p_{14}p_{41} \]
Obtemos a probabilidade de \(A \to \cdots \to A\) usando a regra de adição da Seção 13.4:
\[ P(A \to \cdots \to A) = p_{11}p_{11} + p_{12}p_{21} +p_{13}p_{31} + p_{14}p_{41}=\sum_{i=1}^{4}{p_{1i}p_{i1}} \]
Isso nada mais é do que o produto interno da primeira linha pela primeira coluna da matriz de transição \(P\).
Pelo mesmo argumento, podemos calcular a probabilidade de qualquer outra transição de duas etapas.
O resultado é sempre resumido na fórmula “linha x coluna”. Todas as 16 probabilidades são, portanto, geradas pela multiplicação da matriz \(P^2=PP\).
Podemos verificar se nossa molécula particular de RNA tem a propriedade Markov. Para esse propósito, substituímos as probabilidades de transição desconhecidas pelos valores empíricos e \(P^2\) por \(F^2\). Assim obtemos:
\[ F^2=\left[\begin{array}{rrrr} 0.153 & 0.330 & 0.387 & 0.130\\ 0.149 & 0.335 & 0.361 & 0.155\\ 0.139 & 0.352 & 0.373 & 0.136\\ 0.124 & 0.339 & 0.360 & 0.177 \end{array}\right] \]
Por fim, confrontamos o resultado desse cálculo com as frequências relativas das transições de duas etapas que obtemos por contagem direta da Fig. 14.2:
\[ \left[\begin{array}{rrrr} 0.065 & 0.323 & 0.548 & 0.064\\ 0.153 & 0.444 & 0.208 & 0.195\\ 0.161 & 0.222 & 0.457 & 0.160\\ 0.156 & 0.437 & 0.313 & 0.094 \end{array}\right] \]
A concordância é ruim. Um tratamento estatístico, não reproduzido aqui, revela que existe uma diferença significativa entre as matrizes. Concluímos que nossa molécula de RNA não tem a propriedade de Markov, o que significa que, neste caso, as transições subsequentes dependem umas das outras.
# Matrizes
F2 <- matrix(c(
0.153, 0.330, 0.387, 0.130,
0.149, 0.335, 0.361, 0.155,
0.139, 0.352, 0.373, 0.136,
0.124, 0.339, 0.360, 0.177
), nrow = 4, byrow = TRUE)
O2 <- matrix(c(
0.065, 0.323, 0.548, 0.064,
0.153, 0.444, 0.208, 0.195,
0.161, 0.222, 0.457, 0.160,
0.156, 0.437, 0.313, 0.094
), nrow = 4, byrow = TRUE)
# Frobenius absoluta
F_norm <- sqrt(sum((F2 - O2)^2))
# Normalização
m <- nrow(F2)
n <- ncol(F2)
F_norm_max <- sqrt(m * n)
F_norm_rel <- F_norm / F_norm_max
# Resultado
cat("Distância de Frobenius absoluta:", round(F_norm, 4), "\n")
Distância de Frobenius absoluta: 0.346
Distância relativa normalizada: 0.0865
Interpretação da distância relativa
A teoria da cadeia de Markov desempenha um papel cada vez mais importante em áreas como dinâmica populacional, genética, evolução, ecologia, fisiologia, epidemiologia, comportamento animal.
Um psicólogo realiza um experimento com ratos. Ela coloca um rato em uma gaiola com três compartimentos rotulados como 1, 2 e 3, conforme mostrado na Figura 4.5. Os ratos são treinados para selecionar uma porta aleatoriamente e mover-se de um compartimento para outro quando um sino toca.
Shahin, 2014, p. 293
Questões:
Solução
A matriz de transição \(T\) da cadeia de Markov é construída calculando as probabilidades de mover-se de uma sala para outra.
Shahin, 2014, p. 294
Se o rato está inicialmente na sala 2, o vetor de estado inicial é:
\[ X_0 = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \]
Calculamos a distribuição de estado após três toques do sino multiplicando a matriz de transição \(T\) pelo vetor de estado inicial três vezes:
\[ X_3 = T^3 X_0 \]
A distribuição resultante mostra a probabilidade de estar em cada sala após três toques. Em particular, a probabilidade de estar na sala 1 é aproximadamente 0.3889.
Para encontrar o vetor de estado estacionário \(X_k\), precisamos resolver a equação \(T X = X\) ou o sistema homogêneo \((T - I) X = 0\).
O vetor de estado estacionário é:
\[ X_k \rightarrow \begin{pmatrix} 0.3 \\ 0.4 \\ 0.3 \end{pmatrix} \]
Matriz de Transição T:
Distribuição de Estado Após 3 Passos:
Distribuição de Estado Após 100 Passos:
Vetor de Estado Estacionário:
\(X_k \rightarrow \begin{pmatrix} 0.3 \\ 0.4 \\ 0.3 \end{pmatrix}\)
Indica que, a longo prazo, 30% dos ratos estarão no compartimento 1, 40% no compartimento 2, e 30% no compartimento 3.
Este modelo de Markov fornece uma análise detalhada da movimentação dos ratos entre os compartimentos. Ele permite prever tanto o comportamento de curto prazo (após três toques do sino) quanto o comportamento de longo prazo (vetor estacionário).
# matriz de transição
T <- matrix(c(
0, 1/2, 1/3,
2/3, 0, 2/3,
1/3, 1/2, 0
), nrow = 3, byrow = TRUE)
R <- c("R1", "R2", "R3")
# Nomear as linhas e colunas
rownames(T) <- colnames(T) <- R
# Vetor de estado inicial (começando no compartimento 2)
s0 <- c(0, 1, 0)
# Função para calcular a distribuição de estado após n passos
calcular_estado <- function(T, s0, passos) {
estado_atual <- s0
for (i in 1:passos) {
estado_atual <- T %*% estado_atual
}
return(estado_atual)
}
# Calcular a distribuição de estado após 3 passos
s3 <- calcular_estado(T, s0, 3)
s100 <- calcular_estado(T, s0, 100)
# Resultados
print("Matriz de Transição T:")
[1] "Matriz de Transição T:"
R1 R2 R3
R1 0.0000000 0.5 0.3333333
R2 0.6666667 0.0 0.6666667
R3 0.3333333 0.5 0.0000000
[1] "Distribuição de Estado após 3 Passos:"
[,1]
R1 0.3888889
R2 0.2222222
R3 0.3888889
[1] "Distribuição de Estado após 100 Passos:"
[,1]
R1 0.3
R2 0.4
R3 0.3
# Sistema homogêneo
M <- (T - diag(dim(T)[1]))
sol <- pracma::nullspace(M)
estacionario <- proportions(sol)
rownames(estacionario) <- R
print("Vetor de Estado Estacionário:")
[1] "Vetor de Estado Estacionário:"
[,1]
R1 0.3
R2 0.4
R3 0.3
Proofs involving ordinary least squares: Wikipedia
A álgebra matricial é indispensável em áreas como análise estatística multivariada, planejamento de experimentos e análise de variância e covariância.
Para dar uma ideia de tal aplicação estatística, consideramos um modelo de regressão. Suponha que nos seja dado um diagrama de dispersão consistindo de pontos \((x_i,y_i)\), \(i = 1,2, \ldots, n\), em que \(x_i\) e \(y_i\) são medidas intervalares. Suponha ainda que os pontos estejam próximos de uma certa curva e que queremos ajustar uma função quadrática. a equação é:
\[ y=\beta_0+\beta_1x+\beta_2x^2 \]
em que as constantes \(\beta_0\), \(\beta_1\) e \(\beta_2\) são desconhecidas. As coordenadas \(x_i\) e \(y_i\) dos pontos não satisfazem exatamente uma equação quadrática. Portanto, escrevemos:
\[ Y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\varepsilon_i \]
em que \(\varepsilon_i\) é o termo de erro.
O sistema com \(n\) equações lineares nos parâmetros pode ser escrito na forma matricial:
\[ Y=X\beta+\varepsilon \]
\[ Y=\left[\begin{array}{c} Y_1\\ Y_2\\ \vdots\\ Y_n \end{array}\right] \]
\[ X=\left[\begin{array}{ccc} 1 & x_1 & x_1^2\\ 1 & x_2 & x_2^2\\ \vdots & \vdots & \vdots\\ 1 & x_n & x_n^2 \end{array}\right] \]
\[ \beta=\left[\begin{array}{c} \beta_0\\ \beta_1\\ \beta_2 \end{array}\right] \]
\[ \varepsilon=\left[\begin{array}{c} \varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_n \end{array}\right] \]
A equação matricial \(Y=X\beta+\varepsilon\) é chamada de modelo linear, pois é linear em \(\beta\).
Para encontrar coeficientes adequados, aplicamos o método dos mínimos quadrados (OLS).
Minimizamos a soma dos quadrados dos erros:
\[ \sum_{i=1}^{n}{\varepsilon_i^2}=\varepsilon^{\prime}\varepsilon \]
Isso significa que os coeficientes desconhecidos \(\beta_0\), \(\beta_1\) e \(\beta_2\) são determinados de tal forma que:
\[ Q=\varepsilon^{\prime}\varepsilon=(Y-X\beta)^{\prime}(Y-X\beta) \]
assume o menor valor possível. A solução pode ser encontrada por métodos de álgebra matricial ou podemos aplicar cálculo diferencial.
Conforme Teichroew (1964, p. 580-3) e Tay (2018):
\[ Q=\varepsilon^{\prime}\varepsilon=(Y-X\beta)^{\prime}(Y-X\beta)=\beta^{\prime}X^{\prime}X\beta-2Y^{\prime}X\beta+Y^{\prime}Y \]
Condição de primeira ordem (necessária):
\[ \begin{align} \dfrac{\partial Q}{\partial \beta}&=2X^{\prime}X\beta-2X^{\prime}Y\\ 0&=2X^{\prime}Xb-2X^{\prime}Y\\ b&=(X^{\prime}X)^{-1}X^{\prime}Y \end{align} \]
Condição de segunda ordem (suficiente):
\[ \left.\dfrac{\partial^2 Q}{\partial \beta\partial \beta^{\prime}}\right|_{\beta=b}=2X^{\prime}X \]
A matrix de informação \(XX^{\prime}\) é uma matriz simétrica \(3\times 3\) e também é uma matriz definida positiva (todos os autovalores são positivos), i.e., \(XX^{\prime}>0\). Portanto, \(b\) é a estimativa de mínimos quadrados:
\[ Q_{min}=Y^{\prime}(Y-Xb)=Y^{\prime}e=e^{\prime}e \]
Sendo que \(e=Y-Xb\) é o vetor de resíduos.
Tem-se que:
\[ \begin{align} y&=Xb\\ Y-e&=Xb\\ Y&=Xb+e \end{align} \]
Além disso, a soma dos quadrados dos resíduos é nula:
\[ e^{\prime}1=0 \]
A matriz de planejamento \(X\) e o vetor das estimativas resíduos são ortogonais (\(X\) e \(e\) são matematicamente independentes):
\[ \begin{align} y&=Xb\\ y-Xb&=0\\ X^{\prime}y-X^{\prime}Xb&=0\\ X^{\prime}(y-Xb)&=0\\ X^{\prime}e&=0 \end{align} \]
O estimador OLS de \(\beta\) é não-viesado se \(\mathbb{E}(\varepsilon)=0\):
\[ \begin{align} b&=(X^{\prime}X)^{-1}X^{\prime}Y\\ b&=(X^{\prime}X)^{-1}X^{\prime}(X\beta+\varepsilon)\\ b&=(X^{\prime}X)^{-1}X^{\prime}X\beta+(X^{\prime}X)^{-1}X^{\prime}\varepsilon\\ b&=\beta+(X^{\prime}X)^{-1}X^{\prime}\varepsilon\\ \mathbb{E}(b)&=\beta+(X^{\prime}X)^{-1}X^{\prime}\mathbb{E}(\varepsilon)\\ \mathbb{E}(b)&=\beta \end{align} \]
Os valores preditos de \(Y\) são:
\[ y=Xb=X(X^{\prime}X)^{-1}X^{\prime}Y=PY \]
em que \(P=X(X^{\prime}X)^{-1}X^{\prime}\) é chamada de matriz de projeção (projection matrix) ou matriz chapéu (hat matrix).
\(P\) é simétrica (\(P^{\prime}=P\)) e idempotente (\(P^2=P\)).
O vetor da variável dependente predita \(y\) e o vetor das estimativas dos resídous são ortogonais (\(y\) e \(e\) são matematicamente independentes):
\[y^{\prime}e=0\]
O \(R^2\) é uma medida estatística de quão próximos os pontos observados estão da função de regressão ajustada. Esta medida varia entre 0 (ausência de ajuste) e 1 (ajuste perfeito). \(R^2\) é conhecido como o coeficiente de determinação ou o coeficiente de determinação múltipla para a regressão múltipla.
\[ R^2=\dfrac{\dfrac{y^{\prime}Y}{n}-\bar{Y}^2}{\dfrac{Y^{\prime}Y}{n}-\bar{Y}^2} \]
A estimativa OLS pode ser vista como uma projeção no espaço linear estendido pelos regressores. X1 e X2 referem-se às colunas da matriz de X.) Wikipedia: https://en.wikipedia.org/wiki/Ordinary_least_squares
LeastSquares[{{1, 1}, {1, 2}, {1, 3}}, {5, 8, 25}]
[1] 5 8 25
[,1] [,2]
[1,] 1 1
[2,] 1 2
[3,] 1 3
[,1]
[1,] -7.333333
[2,] 10.000000
[,1]
[1,] -7.333333
[2,] 10.000000
[1] TRUE
[1] TRUE
[,1]
[1,] 32.66667
[,1]
[1,] 32.66667
[,1]
[1,] 0
[,1]
[1,] 0
[,1]
[1,] 0.8595989
plot(X[,2], Y,
main=paste0("OLS: y = Xb\nR^2 = ", round(R2,2)),
ylab="VD",
xlab="X",
ylim=1.1*c(min(y), max(y)))
curve(b[1]+b[2]*x, min(X[,2]), max(X[,2]), add=TRUE)
Y=[5;8;25]
X=[1 1;1 2;1 3]
b=inv(X'*X)*X'*Y
lsq(X,Y)
y=X*b
e=Y-y
P=X*inv(X'*X)*X'
Q=Y'*e
Q=e'*e
esum=e'*[1 1 1]'
ortg=y'*e
n=length(Y)
R2=(y'*Y/n-mean(Y)^2)/(Y'*Y/n-mean(Y)^2)
quit
Scilab 2023.1.0 (May 23 2023, 09:23:00)
Y =
5.
8.
25.
X =
1. 1.
1. 2.
1. 3.
b =
-7.3333333
10.
ans =
-7.3333333
10.000000
y =
2.6666667
12.666667
22.666667
e =
2.3333333
-4.6666667
2.3333333
P =
0.8333333 0.3333333 -0.1666667
0.3333333 0.3333333 0.3333333
-0.1666667 0.3333333 0.8333333
Q =
32.666667
Q =
32.666667
esum =
-1.066D-14
ortg =
-1.492D-13
n =
3.
R2 =
0.8595989
Em cada repetição, avalia-se se o modelo de regressão linear polinomial interpola exatamente os \(n\) pontos (isto é, se os valores ajustados coincidem com os valores observados de \(y\), até o máximo erro numérico absoluto menor que 1e-5).
# Gerar n_pts pontos com x distintos
n_pts <- 4
x <- sort(runif(n_pts, min = 1, max = 10))
y <- runif(n_pts, min = -5, max = 5)
# Ajustar polinômio de grau n - 1
modelo <- lm(y ~ poly(x, degree = n_pts-1, raw = TRUE))
# Predição nos pontos de entrada
pred <- predict(modelo, newdata = data.frame(x = x))
tol <- 1e-5
# Verificar se o ajuste é exato (erro numérico próximo de zero)
erro <- max(abs(pred - y)) < tol
# Plotar
plot(x, y, pch = 16, col = "black", main = "Polinômio de grau n-1 ajustado a n pontos", ylim=2*c(min(y),max(y)))
curve(predict(modelo, newdata = data.frame(x = x)),
from = min(x), to = max(x), col = "blue", lwd = 2, add = TRUE)
abline(h=0, lty=2)
Ajuste exato? TRUE
Definimos vetores como matrizes linha e coluna. Vetores com dois ou três elementos podem ser interpretados geometricamente. Por um lado, esta interpretação permite uma melhor compreensão das operações matriciais. Por outro lado, vetores interpretados geometricamente podem ser aplicados a vários problemas das ciências da vida.
Primeiro consideramos um sistema retangular de coordenadas \(xy\) no plano.
Seja \((x, y)\) um ponto arbitrário no plano. Ao ponto \((x, y)\) associamos um segmento de reta direcionado, também chamado de seta.
Sua cauda está na origem \(O\) do sistema de coordenadas e a ponta coincide com o ponto \((x, y)\). Seja a o vetor coluna
\[ a=\left[\begin{array}{c} x\\ y \end{array}\right] \]
Então há uma correspondência biunívoca entre o vetor \(a\) e a seta descrita acima. É comum usar a palavra “vetor” tanto para a matriz \(a\) quanto para a seta (ver Fig. 14.3). Chamamos \(x\) e \(y\) de coordenadas ou componentes do vetor \(a\).
Fig. 14.3. Um vetor a é representado por uma seta com coordenadas x, y.
a <- c(2,1.5)
xlim <- c(0,2.5)
ylim <- c(0,2)
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1)
abline(v=0, h=0, col="gray")
matlib::vectors(a)
O significado original de “vetor” era uma quantidade direcionada, como uma força, uma velocidade ou uma aceleração. As quantidades não direcionadas são chamadas de escalares. Por exemplo, massa, densidade e temperatura são escalares. Mais tarde, a palavra “vetor” também foi usada no sentido de um segmento de linha direcionado ou uma seta. Nas últimas décadas, o significado mudou novamente, ou seja, para matrizes de linha e coluna. Da mesma forma, podemos chamar \((x, y)\) uma matriz, um vetor ou um ponto. Na álgebra matricial, os vetores coluna são preferidos aos vetores linha.
Sejam dois vetores \(a_1\) e \(a_2\):
Fig. 14.4. A soma de dois vetores a1 e a2 é um vetor que pode ser construído por meio de um triângulo. Os componentes de a1 + a2 são x1 + x2 e y1 + y2.
[,1]
[1,] 48.94519
xlim <- c(0,6)
ylim <- c(0,5)
# proper geometry requires asp=1
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1,
main = expression(theta == 49))
abline(v=0, h=0, col="gray")
matlib::vectors(rbind(a,b), col=c("red", "blue"), cex.lab=c(2, 2))
text(.5, .37, expression(theta), cex=1.5)
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1,
main = expression(theta == 49))
abline(v=0, h=0, col="gray")
sum <- a+b
sum
[1] 4.5 4.0
matlib::vectors(rbind(a,b, "a+b"=sum),
col=c("red", "blue", "black"),
cex.lab=c(2, 2, 2.2))
matlib::vectors(sum, origin=a, col="red", lty=2)
matlib::vectors(sum, origin=b, col="blue", lty=2)
Agora tentamos encontrar um significado geométrico para a soma vetorial:
\[ a_1+a_2=\left[\begin{array}{c} x_1+x_2\\ y_1+y_2 \end{array}\right] \]
Na Fig. 14.4, os três vetores \(a_1\), \(a_2\) e \(a_1 + a_2\) são representados. Também é mostrado como um \(a_1+ a_2\) pode ser construído. Para isso, movemos a flecha \(a_2\) temporariamente, paralela a si mesma, de modo que sua cauda coincida com a ponta de \(a_1\).
Agora as duas flechas formam uma linha tracejada.
A soma \(a_1 +a_2\) é representada pela seta que une a origem \(O\) com a ponta da seta deslocada \(a_2\). Diz-se que os três vetores formam um triângulo vetorial. Devido à lei comutativa da adição de matrizes, também poderíamos mover a seta \(a_1\) de modo que ela forme uma linha tracejada com a seta \(a_2\).
É fácil generalizar a adição de vetores para três ou mais vetores.
Se tivermos que adicionar \(a_1, a_2, \ldots , a_n\), formamos uma linha quebrada começando com \(a_1\), digamos, então continuando com as setas deslocadas \(a_2, a_3, \ldots, a_n\). A seta que une a origem \(O\) com a ponta de \(a_n\) representa o vetor \(a_1 + a_2 + \ldots + a_n\). Em aplicações, a soma é frequentemente chamada de vetor resultante e os vetores únicos \(a_1, a_2, \ldots , a_n\) e seus componentes.
A figura que consiste em \(a_1, a_2, \ldots , a_n\) e o vetor resultante é chamada de polígono vetorial. É uma generalização do triângulo vetorial.
Às vezes é preferível construir a soma de dois vetores por meio de um paralelogramo. A Fig. 14.5 explica o procedimento.
Fig. 14.5. A construção de uma soma vetorial por meio de um paralelogramo. A linha P1P2 é paralela à linha OP2 e a linha P2P3 é paralela à linha OP1.
No entanto, esta construção não é prática para mais de dois vetores.
Agora consideramos um múltiplo de um vetor a. O significado geométrico de \(a+a=2a\), \(a+a+a=3a\) etc. é fácil de entender (Fig. 14.6). Também podemos multiplicar um vetor por um número fracionário, digamos por 2/3, ou por um número negativo, digamos por -2.
Obtemos:
\[ ka=k\left[\begin{array}{c} x\\ y \end{array}\right]=\left[\begin{array}{c} kx\\ ky \end{array}\right] \]
Ambas as coordenadas, \(x\) e \(y\), devem ser multiplicadas por \(k\) (Fig. 14.6).
Fig. 14.6. Multiplicação de um vetor II por um número k. Se multiplicamos um vetor por um número positivo, sua direção permanece inalterada. No entanto, a multiplicação por um número negativo produz um vetor apontando na direção oposta.
Ambas as coordenadas, \(x\) e \(y\), devem ser multiplicadas por \(k\) (Fig. 14.6).
A subtração é definida por
\[ a - b= a +(- b) \]
ou seja, a subtração de \(b\) é equivalente à adição do vetor oposto \(-b\).
A Fig. 14.7 mostra a construção e o significado geométrico de \(a-b\).
Fig. 14.7. Os três vetores a, b e a-b formam um triângulo vetorial tal que b + (a - b) = a.
Notar que
\[ a-a=\left[\begin{array}{c} x\\ y \end{array}\right]-\left[\begin{array}{c} x\\ y \end{array}\right]= \left[\begin{array}{c} 0\\ 0 \end{array}\right] \]
O resultado é chamado de vetor zero. Tanto a ponta quanto a cauda desse vetor caem na origem \(O\). Esse vetor em particular não tem direção.
Observe que não precisamos fazer uma suposição sobre as unidades nas quais \(x\) e \(y\) são medidos. Todos os resultados anteriores permanecem válidos se, por exemplo, \(x\) for medido em gramas e \(y\) em segundos. Isso não será mais verdadeiro quando introduzirmos o valor absoluto de um vetor. De agora em diante, assumimos que ambas as coordenadas \(x\) e \(y\) são medidas na mesma unidade.
Definimos o valor absoluto (norma-\(L^2\), magnitude) de um vetor \(a^{\prime} = [x\; y]\) por
\[ |a|=\sqrt{a^{\prime}a}=\sqrt{x^2+y^2} \]
Sendo \(a^{\prime}a\) chamado de produto interno.
Se \(a=0\), \(|a|=0\); caso contrário, o valor absoluto é um número positivo.
[,1]
[1,] 45
xlim <- c(0,1.5)
ylim <- c(0,1.5)
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1,
main = expression(theta == 45))
abline(v=0, h=0, col="gray")
matlib::vectors(rbind(a,b), col=c("red", "blue"), cex.lab=c(2, 2))
text(.25, .10, expression(theta), cex=1.5)
[1] 1.414214
[1] 1
[1] 1
No caso especial em que \(x\) e \(y\) são medidos em uma unidade de comprimento, então \(|a|\) é a distância de \(O\) ao ponto \((x, y)\) em virtude do teorema de Pitágoras (Fig. 14.3). Nesse caso, \(|a|\) também é chamado de comprimento do vetor \(a\).
Fig. 14.8. Prova da fórmula pitagórica.
A lei comutativa funciona para o produto interno. Como já sabemos a lei comutativa em geral não é válida para a multiplicação de matrizes.
Para abordar a geometria do produto interno, consideramos primeiro um caso especial.
Assumimos que as setas \(a_1\) e \(a_2\) são perpendiculares (Fig. 14.8). Os vetores \(a_1\), \(a_2\) e \(a_1-a_2\) formam um triângulo retângulo, i.e., são perpendiculares ou ortogonais.
Segue do teorema de Pitágoras que:
\[ |a_1-a_2|^2=|a_1|^2+|a_2|^ \]
Note que, de modo geral (Fig. 14.9):
\[ |a_1-a_2|^2=|a_1|^2+|a_2|^2-2a_{1}^{\prime}a_2 \]
Se \(a_1\) e \(a_1\) são perpendiculares, então \(a_{1}^{\prime}a_2=0\).
Fig. 14.9. Prova da fórmula não pitagórica.
Esta equação vale se, e somente se, os dois vetores diferentes de zero \(a_1\) e \(a_2\) são perpendiculares ou ortogonais. A condição \(a_{1}^{\prime}a_2=0\) é, portanto, chamada de condição de ortogonalidade.
O resultado pode ser estendido para o caso em que \(a_1\) e \(a_2\) formam um ângulo arbitrário \(\alpha\) (Fig. 14.9). Temos que substituir o teorema de Pitágoras pela lei dos cossenos que afirma em nossa notação:
\[ |a_1-a_2|^2=|a_1|^2+|a_2|^2-2|a_{1}||a_2|\cos (\alpha) \]
Portanto, o produto escalar (dot product) dos dois vetores é expresso por:
\[ a_{1}^{\prime}a_2=|a_{1}||a_2|\cos (\alpha) \]
Produto escalar de vetores. Percebe-se que |A|cos(θ) é a projeção escalar de A em B.
Em palavras: O produto interno de dois vetores \(a_1\) e \(a_2\) é igual ao produto de três fatores: os valores absolutos de \(a_1\) e \(a_2\) e o cosseno do ângulo entre os dois vetores.
u <- c(3,1)
v <- c(1,3)
xlim <- c(0,4)
ylim <- c(0,4)
plot(xlim, ylim, type="n", xlab="x", ylab="y", asp=1)
abline(v=0, h=0, col="gray")
matlib::vectors(rbind(u,v),
col=c("red", "blue"),
cex.lab=c(2, 2))
# projection of vectors
O <- c(0,0)
P <- matlib::Proj(v,u)
matlib::vectors(P, labels="P", lwd=3)
matlib::vectors(v, origin=P)
matlib::corner(O, P, v)
[1] 3.162278
[1] 1.897367
Se um dos vetores for o vetor zero, o produto é zero, e \(\alpha\) é indeterminado. Se o ângulo \(\alpha\) é agudo, então \(\cos(\alpha)> 0\), e o produto (14.4.13) é positivo. Se, no entanto , o ângulo \(\alpha\) é obtuso, então \(\cos(\alpha)< 0\), e o produto (14.4.13) é negativo. Já discutimos o caso \(\alpha=90^{\circ}\) que implica \(\cos(\alpha)= 0\).
A fórmula \(a_{1}^{\prime}a_2=|a_{1}||a_2|\cos (\alpha)\) é frequentemente usada para calcular o ângulo \(\alpha\) entre dois vetores diferentes de zero. Obtemos:
\[ \cos(\alpha)=\dfrac{a^{\prime}_{1}a_2}{|a_1||a_2|} \]
De fato, a extensão dos resultados do espaço bidimensional para o tridimensional é tão simples que os matemáticos não puderam deixar de inventar o espaço quadridimensional (Fig. 14.10). A falta de intuição geométrica não era barreira. A passagem do espaço tridimensional para o quadridimensional era tão tentadora que precisava ser empreendida. Convidamos o leitor a realizar esta etapa sem receber nenhuma ajuda adicional. (Veja também o Exemplo 14.5.7).
Fig. 14.10. Um vetor tridimensional a com coordenadas x, y, z.
vec <- rbind(diag(3), c(1,1,1))
rownames(vec) <- c("x", "y", "z", "a")
invisible(rgl::open3d())
matlib::vectors3d(vec, color=c(rep("black",3), "black"), lwd=2)
# draw the XZ plane, whose equation is Y=0
rgl::planes3d(0, 0, 1, 0, col="gray", alpha=0.2)
matlib::vectors3d(c(1,1,0), col="black", lwd=2)
# show projections of the unit vector J
rgl::segments3d(rbind(c(1,1,1), c(1, 1, 0)))
rgl::segments3d(rbind(c(0,0,0), c(1, 1, 0)))
rgl::segments3d(rbind(c(1,0,0), c(1, 1, 0)))
rgl::segments3d(rbind(c(0,1,0), c(1, 1, 0)))
# show some orthogonal vectors
p1 <- c(0,0,0)
p2 <- c(1,1,0)
p3 <- c(1,1,1)
p4 <- c(1,0,0)
matlib::corner(p1, p2, p3, col="black")
#matlib::corner(p1, p4, p2, col="blue")
#matlib::corner(p1, p4, p3, col="black")
rgl::rglwidget()
Não só foi possível trabalhar com vetores no espaço quadridimensional com facilidade, como as aplicações nas ciências naturais também se tornaram bastante bem-sucedidas. A generalização posterior da geometria de vetores para espaços \(n\)-dimensionais com \(n = 5, 6, \ldots\) era uma questão de rotina. Por volta de 1900, até espaços infinitamente dimensionais foram inventados.
Hoje, a álgebra vetorial em espaços \(n\)-dimensionais é uma ferramenta útil na análise estatística multivariada, especialmente na análise de variância e covariância.
Na Fig. 14.11, um corpo humano está em repouso sobre um plano inclinado. Qual é a força que está tentando puxar o corpo para baixo ao longo do plano e qual é a força que pressiona o corpo em direção ao plano?
Seja \(F_1\) a força que puxa o corpo para baixo ao longo do plano e \(F_2\) a força perpendicular ao plano que pressiona o corpo em direção ao plano. Em mecânica, aprendemos que essas duas forças são causadas pela gravitação e que a soma de seus vetores deve ser a força gravitacional total denotada por \(F\). Portanto:
\[ F_1+F_2=F \]
Fig. 14.11. Um corpo em um plano inclinado. A adição vetorial relaciona as forças F1 e F1 com a força F da gravitação.
As forças \(F_1\) e \(F_2\) são chamadas de componentes de \(F\). Seja \(\alpha\) o ângulo de inclinação do plano. Como \(F_2\) é perpendicular ao plano, \(\alpha\) também é o ângulo entre \(F_2\) e \(F\). Da Fig. 14.11 obtemos as seguintes fórmulas para as magnitudes dos componentes de \(F\):
\[ |F_1|=|F|\sin (\alpha)\\ |F_2|=|F|\cos (\alpha) \]
Por exemplo, se \(\alpha=30^{\circ}\), obtemos \(\sin (\alpha)=0.500\), \(\cos (\alpha)=0.866\). Assim \(F_1\) é apenas 50% e \(F_2\) apenas 86.6% do peso corporal. Por um lado, para impedir o deslizamento do corpo, uma força de magnitude \(\frac{1}{2}F_1\) tem que agir na direção oposta a \(F_1\). Por outro lado, a pressão contra o plano é aliviada em 13.4% (=100% - 86.6%) em comparação com o corpo deitado um plano horizontal.
Como exemplo de alavanca, tomamos o antebraço (Fig. 14.13). Seu ponto de apoio é o cotovelo. Quando o ângulo entre o braço e o antebraço não é de 90°, a força \(F\) gerada pelo flexor do antebraço deve ser decomposta em um componente \(F_1\) perpendicular ao antebraço e um componente \(F_2\) paralelo ao antebraço. A força \(F_2\) não gera nenhuma rotação do antebraço. Apenas \(F_1\) atua na direção correta. Se o ângulo entre o braço e o antebraço se aproxima de 180°, a magnitude de \(F_1\) é consideravelmente menor que a de \(F\). Assim, parte da força muscular é “perdida”. Por álgebra vetorial podemos escrever:
\[ F_1=F-F_2 \]
Fig. 14.12. Forças atuando em um osso quebrado. A força F pode ser causada pela gravidade ou por um músculo. Chamamos F1 de força de cisalhamento.
Fig. 14.13. A força F1 que levanta o antebraço pode ser consideravelmente menor em magnitude do que a força F gerada pelo flexor do antebraço.
A Fig. 14.14 mostra como uma perna pode ser esticada por uma polia para fins terapêuticos. Denotamos por \(F_1\) a força vertical do peso. O fio da polia tem a mesma tensão em todos os lugares. Portanto, as forças \(F_2\) e \(F_3\) na Fig. 14.14a ou b têm a mesma magnitude que \(F_1\), ou seja,
\[ |F_1|=|F_2|=|F_3| \]
A força de alongamento \(F\) é a resultante de \(F_2\) e \(F_3\). Então:
\[ F=F_2+F_3 \]
Quando o ângulo entre \(F_2\) e \(F_3\) tende a zero, então \(|F|\) aumenta e tende a \(|F_2|+|F_3|=2|F_1|\) como é facilmente visto na Fig. 14.14a. Quando, entretanto, o ângulo entre \(F_2\) e \(F_3\) tende a 180°, \(|F|\) diminui e tende a zero (Fig. 14.14b). Então:
\[ 0 \le |F| \le 2|F_1| \]
Fig. 14.14. Tração da perna com peso. A magnitude da força resultante depende do ângulo entre as forças F2 e F3.
Um corpo consiste em um número muito grande de partículas que estão sujeitas à gravidade. Para encontrar o centro de gravidade, começamos com o caso simples de dois pontos \(P_1\) e \(P_2\) de massas iguais (Fig. 14.15).
Então, por simetria, o centro de gravidade é o ponto médio de \(P_1\) e \(P_2\). Denotamos por \(C\). Para encontrar o ponto médio \(C\) por álgebra vetorial, denotamos as setas com cauda na origem \(O\) e pontas em \(P_1\) e \(P_2\) por \(a_1\) e \(a_2\) respectivamente. Usando o paralelogramo da Fig. 14.5, determinamos a soma vetorial de \(a_1\) e \(a_2\). Como as diagonais de um paralelogramo se dividem, o vetor \(c\) que aponta para \(C\) é:
\[ c=\dfrac{a_1+a_2}{2} \]
Fig. 14.15. O ponto médio de dois pontos definidos pelos vetores a1 e a2.
Isso é simplesmente a média aritmética dos dois vetores \(a_1\) e \(a_2\).
Se \(a_{1}^{\prime}=[x_1\; x_2]\) e \(a_{2}^{\prime}=[y_1\; y_2]\), então:
\[ c^{\prime}=\left[\dfrac{x_1+x_2}{2} \; \dfrac{y_1+y_2}{2}\right] \]
Se nos forem dados \(n\) pontos de mesma massa pelos vetores \(a_1, a_2, \ldots, a_n\), a mecânica mostra que o centro de gravidade é dado pelo vetor
\[ c=\dfrac{\sum_{i=1}^{n}{a_i}}{n} \]
A mesma fórmula é de importância básica na análise estatística das direções aplicadas, por exemplo, a pássaros migratórios e homing. Para obter detalhes, consulte Batschelet (1965).
Uma generalização adicional ocorre quando diferentes massas \(M_i\) estão localizadas nos pontos dados pelos vetores \(a_i\) (\(i = 1,2, \ldots , n\)). Nesse caso, a média aritmética ordinária se transforma na média aritmética ponderada
\[ c=\dfrac{\sum_{i=1}^{n}{M_i a_i}}{\sum_{i=1}^{n}{M_i}}=\sum_{i=1}^{n}{\dfrac{M_i }{\sum_{i=1}^{n}{M_i}}}a_i \]
A Fig. 14.16 ilustra o caso de três massas. A fórmula é válida no espaço bidimensional e tridimensional. Também é útil em espaços de dimensões superiores, embora perca seu significado físico. Uma aplicação importante de (14.5.9) é feita na análise estatística multivariada.
Fig. 14.16. O centro de gravidade C de três massas M1, M2, M3 localizadas nos pontos dados pelos vetores a1, a2, a3, respectivamente.
Newton's laws of motion
: Wikipedia
Existem inúmeras aplicações biológicas de vetores no estudo da locomoção animal. Vamos nos limitar a discutir o salto dos gafanhotos (Fig. 14.15).
Fig. 14.17. gafanhoto pulando.
Chamamos o intervalo de tempo desde o início da ação muscular até que os pés deixem o solo de decolagem ou período de aceleração.
Durante este período, o centro de gravidade move-se de um ponto \(C\) para um certo ponto \(C^{\prime}\).
Os dois pontos \(C\) e \(C^{\prime}\) formam um segmento de reta direcionado ou um vetor que denotamos por \(s\).
Seja \(\alpha\) o ângulo entre \(s\) e a linha horizontal. Então \(|s| \sin (\alpha)\) é a componente vertical de \(s\).
Durante o período de decolagem, a força que acelera o corpo do gafanhoto não permanece constante em magnitude, mas podemos operar com uma força média que denotamos por \(F\). Essa força \(F\) é chamada de empuxo.
Agora queremos calcular a altura máxima \(h\) sobre o solo atingida pelo gafanhoto saltador. Por um lado, a parte da energia cinética que é transformada em energia potencial é a componente vertical de \(F\) vezes a componente vertical de \(s\) é expressa por:
\[ W_c=|F|\sin (\alpha)|s|\sin (\alpha)=|F||s|\sin^2 (\alpha) \]
Sendo que \(W_c\) é o trabalho cinético realizado por uma força \(F\) sobre um objeto ao longo de uma distância \(s\) paralela a \(F\).
Por outro lado, o trabalho potencial alcançado no ponto mais alto é
\[ W_p=m\times g\times h \]
Sendo que \(W_p\) é o trabalho (potencial) realizado por uma força \(F\) sobre um objeto ao longo de uma distância \(h\) paralela a \(f=m \times g\), \(m\) é a massa do animal (em gramas) e \(g\) a aceleração da gravidade (= 9.81 m seg-2). Lembra que a segunda lei de Newton estabelece que \(f=m \times a\), em que \(a\) é aceleração.
Se negligenciarmos a pequena perda de energia devido à resistência do ar, igualamos as duas energias. Então dividindo por \(m \times g\) nós obtemos
\[ \begin{align} W_c&=W_p\\ m\times g\times h&=|F||s|\sin^2 (\alpha) \end{align} \]
Esta fórmula também pode ser usada para calcular o \(|F|\) quando as outras grandezas são conhecidas a partir das medições.
\[ |F| = \dfrac{m\times g\times h}{|s|\sin^2 (\alpha)} \]
Frequências gênicas no grupo sanguíneo ABO têm sido investigadas em várias populações. Se distinguirmos entre quatro alelos \(A_1, A_2 , B, O\), as seguintes frequências relativas, \(f_{ki}\) foram relatadas (Cavalli-Sforza e Edwards, 1967, p. 250):
\[ \begin{array}{c|c|c|c|c} \text{Alelo} & \text{Esquimó} & \text{Bantu} & \text{Inglês} & \text{Coreano} \\ & f_{1i} & f_{2i} & f_{3i} & f_{4i} \\ \hline A_1 & 0.2914 & 0.1034 & 0.2090 & 0.2208 \\ A_2 & 0.0000 & 0.0866 & 0.0696 & 0.0000 \\ B & 0.0316 & 0.1200 & 0.0612 & 0.2069 \\ O & 0.6770 & 0.6900 & 0.6602 & 0.5723 \\ \hline \text{Total} & 1.0000 & 1.0000 & 1.0000 & 1.0000 \\ \end{array} \]
Surge a pergunta:
Quão perto está uma população de outra?
Em outras palavras: devemos encontrar uma medida adequada para uma distância genética.
Um método usando álgebra vetorial foi proposto por Cavalli-Sforza e Edwards (1967). Num primeiro passo representamos cada população por um vetor unitário, ou seja, um vetor de valor absoluto 1. Para isso tomamos a raiz quadrada de cada frequência, digamos
\[ x_{ki}=\sqrt{f_{ki}} \]
Se \(\sum_{i=1}^{4}{f_{ki}}=1\), então para cada uma das quatro populações, tem-se que \(\sum_{i=1}^{4}{x^{2}_{ki}}=1\).
Isto significa que cada um dos quatro vetores é um vetor unitário, i.e., \(|a_k|=1\), sendo que \(a_k^{\prime}=[x_1\; x_2\; x_3\; x_4]\), \(i=1\; (\text{esquimó}), 2\; (\text{bantus}), 3\; (\text{inglês}), 4\; (\text{coreano})\).
Os bantus ou bantos constituem um grupo etnolinguístico localizado principalmente na África subsariana e que engloba cerca de 400 subgrupos étnicos diferentes.
No espaço quadridimensional, as pontas desses vetores estão localizadas em uma hiper-esfera de raio 1.
Agora é plausível usar o ângulo entre dois vetores como medida da distância entre as populações correspondentes. Se denotarmos o ângulo entre \(a_1\) (esquimó) e \(a_2\) (bantus), então:
\[ \cos(\theta) = a_{1}^{\prime}a_{2} \]
Note que \(|a_1|=|a_2|=1\).
Aqui está um exemplo numérico:
\[ a_1^{\prime}=[\sqrt{0.2914}\; \sqrt{0.0000}\; \sqrt{0.0316}\; \sqrt{0.6770}]\\ a_2^{\prime}=[\sqrt{0.1034}\; \sqrt{0.0866}\; \sqrt{0.1200}\; \sqrt{0.6900}] \]
\[ a_1^{\prime}=[0.5398\; 0.0000\; 0.1778\; 0.8228]\\ a_2^{\prime}=[0.3216\; 0.2943\; 0.3464\; 0.8307] \]
\[ \begin{align} \cos(\theta) &= a_{1}^{\prime}a_{2}=0.9187\\ \theta&=23.2^{\circ} \end{align} \]
Da mesma forma, encontramos as seguintes distâncias genéticas:
\[ \begin{array}{c|c|c|c} & \text{Bantus} & \text{Inglês} & \text{Coreano} \\ \hline \text{Esquimó} & 23.2^\circ & 16.4^\circ & 16.8^\circ \\ \text{Bantus} & & 9.8^\circ & 20.4^\circ \\ \text{Inglês} & & & 19.6^\circ \end{array} \]
A menor distância genética entre diferentes populações ocorre em bantus-inglês, o maior em esquimó-bantus.
# Criar matriz de distâncias
dist_matrix <- matrix(c(
0, 23.2, 16.4, 16.8,
23.2, 0, 9.8, 20.4,
16.4, 9.8, 0, 19.6,
16.8, 20.4, 19.6, 0),
nrow = 4, byrow = TRUE)
# Nomear linhas e colunas
dimnames(dist_matrix) <- list(
c("Esquimó", "Bantus", "Inglês", "Coreano"),
c("Esquimó", "Bantus", "Inglês", "Coreano")
)
# Converter para objeto de distâncias
d <- as.dist(dist_matrix)
# Cluster hierárquico
hc <- hclust(d)
# Dendrograma
plot(hc, main = "Cluster de Grupos Linguísticos", xlab = "", ylab = "Distância", sub = "")
Cavalli-Sforza & Edwards, 1967, p. 254
Considere uma superfície representando uma função escalar de duas variáveis, como a elevação do terreno dada por \(f(x, y)\), onde \(x\) e \(y\) são coordenadas no plano e \(f(x, y)\) é a altura nesse ponto.
O vetor gradiente \(\nabla f(x, y)\) em um ponto dessa superfície tem:
Interpretação gráfica:
Imagine que você está parado em uma montanha. O vetor gradiente naquele ponto:
Se você seguir o vetor gradiente, estará fazendo o caminho de subida mais rápida possível.
Exemplo:
Para \(f(x, y) = x^2 + y^2\), o gradiente é
\[ \nabla f(x, y) = (2x, 2y) \]
Esse vetor cresce radialmente a partir da origem. A superfície é um paraboloide, e o gradiente aponta para fora, com magnitude proporcional à distância até o centro.
Se quiser, posso gerar um gráfico em R mostrando as curvas de nível e os vetores gradiente sobrepostos. Deseja?
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL","pt_BR.UTF-8"))
# Função escalar
f <- function(x, y) x^2 + y^2
# Gradiente: ∇f = (2x, 2y)
grad_f <- function(x, y) c(2 * x, 2 * y)
# Grade para gráfico 2D
x_seq <- seq(-2, 2, length.out = 10)
y_seq <- seq(-2, 2, length.out = 10)
grid <- base::expand.grid(x = x_seq, y = y_seq)
grid$z <- base::with(grid, f(x, y))
# Vetores gradiente
gradientes <- base::t(base::mapply(grad_f, grid$x, grid$y))
grid$u <- gradientes[, 1]
grid$v <- gradientes[, 2]
# Gráfico 2D com campo vetorial
g <- ggplot2::ggplot(grid, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z), alpha = 0.4) +
ggplot2::geom_segment(ggplot2::aes(xend = x + u * 0.2,
yend = y + v * 0.2),
arrow = ggplot2::arrow(length = grid::unit(0.15, "cm")),
color = "blue") +
ggplot2::scale_fill_gradient(low = "white", high = "black") +
ggplot2::labs(title = "Vetor Gradiente de f(x, y) = x² + y²",
x = "x", y = "y", fill = "f(x, y)") +
ggplot2::coord_fixed() +
ggplot2::theme_minimal()
print(g)
# Grade para gráfico 3D
x_vals <- seq(-2, 2, length.out = 100)
y_vals <- seq(-2, 2, length.out = 100)
grid_3d <- base::expand.grid(x = x_vals, y = y_vals)
grid_3d$z <- base::with(grid_3d, f(x, y))
# Matriz z para superfície
z_matrix <- base::matrix(grid_3d$z, nrow = 100)
# Gráfico 3D interativo com plotly
fig <- plotly::plot_ly(x = ~x_vals, y = ~y_vals, z = ~z_matrix)
fig <- plotly::add_surface(fig,
contours = list(
z = list(
show = TRUE,
usecolormap = TRUE,
highlightcolor = "#ff0000",
project = list(z = TRUE)
)
))
fig <- plotly::layout(fig,
scene = list(
xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "f(x, y)"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = 0.64)
)
),
title = "Superfície 3D da Função f(x, y) = x² + y²")
print(fig)
Álgebra linear é o ramo da matemática que estuda vetores, espaços vetoriais, matrizes e transformações lineares. Seu foco está em estruturas lineares e nas relações entre elas.
Conceitos fundamentais:
Aplicações:
Álgebra linear fornece a linguagem e as ferramentas fundamentais para modelar e resolver problemas em quase todas as áreas quantitativas.
O primeiro passo é definir a anatomia da matriz que será objeto de análise. A matriz quadrada, i.e., a matriz com o número de linhas igual ao número de colunas constitui a partir desse momento o objeto de análise.
Matematicamente, tem-se que a matriz quadrada real \(\underset{k\times k}{A}\), sendo \(k=2,3,\ldots\), consiste em uma tabela cujos \(k^2\) valores numéricos são valores reais, i.e., \(a_{i,j}\in\mathbb{R}\), estão indexados por um número de linha \(i=1,2,\ldots\) e por um número de coluna \(j=1,2,\ldots\).
A diagonal principal de \(\underset{k\times k}{A}\) consiste nos valores de índice \(i=j\), sendo \(i=1,2,\ldots\) e é expressa por \(\text{diag}[a_{1,1},a_{2,2}, \ldots, a_{k,k}]\).
A matrix transposta de matriz diagonal é igual à própria matriz diagonal.
Transpose[DiagonalMatrix[{a,b}]]
Uma matriz quadrada real \(\underset{k\times k}{A}\) é não singular se seu determinante é não nulo, i.e., se \(\det(A)\ne 0\) ou, equivalentemente, se seu posto (rank) é completo, i.e., se \(\text{rank}(A)=k\). Isso garante que \(A\) é invertível.
Se a matriz quadrada tem inversa, pode ser usada em sistema linear \(Ax = b\).
A simétrica tem propriedades especiais (autovalores reais, diagonalizável por matriz ortogonal).
[,1] [,2]
[1,] 1 2
[2,] 3 4
[1] TRUE
[1] FALSE
[1] TRUE
[1] -2
[1] 2
[,1] [,2]
[1,] 2 1
[2,] 1 3
[1] TRUE
[1] TRUE
[1] TRUE
[1] 5
[1] 2
[,1] [,2]
[1,] 1 2
[2,] 2 4
[1] TRUE
[1] TRUE
[1] FALSE
[1] 0
[1] 1
Scilab 2023.1.0 (May 23 2023, 09:23:00)
ans =
-2.
ans =
2.
ans =
5.
ans =
2.
ans =
0.
ans =
1.
Sejam \(A\) uma matriz quadrada de números reais e \(I\) uma matriz identidade quadrada de mesma dimensão de \(A\).
O autovalor, também conhecido por valor próprio, número próprio, raiz característica, é obtido pela solução da seguinte equação:
\[\det(A-\lambda I)=0\]
Pelo teorema espectral da álgebra linear, os autovalores \(\lambda_i\) são números reais se a matriz quadrada \(A\) é simétrica.
A equação anterior gera um polinômio de ordem \(k\) e é denominado polinômio característico.
Os \(k\) autovalores ou raízes características extraídas do polinômio característico formam o espectro ordenado da matriz \(A\), isto é, \(\text{spec}(A)=(\lambda_1, \lambda_2, \ldots, \lambda_k)\), sendo \(\lambda_1\ge \lambda_2\ge \cdots \ge \lambda_k\).
CharacteristicPolynomial[{{a,b}, {c,d}}, lambda] = 0
a^2 + 4 b c - 2 a d + d^2 = 4 b c + (a - d)^2
CharacteristicPolynomial[{{a,b}, {b,d}}, lambda] = 0
a^2 + 4 b^2 - 2 a d + d^2 = (a - d)^2 + 4 b^2 > 0
Uma matriz quadrada \(A\) é simétrica se \(a_{i,j}=a_{j,i}\), i.e., \(A^{\prime}=A\).
Se a matriz \(A\) é simétrica, então seu espectro ordenado é real.
Valores singulares de \(A\) são definidos como:
\[ \sigma_i = \sqrt{\lambda_i(A^{\prime} A)} \]
ou seja, são as raízes quadradas dos autovalores da matriz simétrica positiva semi-definida \(A^{\prime} A\). Por isso, são sempre reais e não-negativos.
Os autovalores podem ser negativos, nulos ou complexos. Os valores singulares sempre são reais não negativos, refletindo a “intensidade” com que a matriz transforma vetores, e são mais estáveis numericamente (robustos).
Exemplo:
[1] 5.3722813 -0.3722813
[1] 5.4649857 0.3659662
[1] 5.4649857 0.3659662
SingularValueList[{{a,b}, {b,d}}]
A matriz inversa da matriz quadrada real \(A\) é denotada por \(A^{-1}\) tal que \(AA^{-1}=I\).
Se a matriz \(A\) é quadrada e não singular, isto é, se \(\text{rank}(A)=k\), então existe e é única a matriz inversa \(A^{-1}\).
A matriz inversa nesse caso é dado por \(A^{-1}=\dfrac{\text{cof}^{\prime}(A)}{\det(A)}\), sendo que \(\text{cof}^{\prime}(A)\) é a matriz de cofatores transposta. Se a matriz \(A\) é simétrica, então \(A^{-1}=\dfrac{\text{cof}(A)}{\det(A)}\).
Se a matriz \(A\) é diagonal, então sua inversa é dada pelo inverso de cada elemento da diagonal.
Cofactor[{{a,b}, {c,d}}]/Det[{{a,b}, {c,d}}]
Inverse[Cofactor[{{a,b}, {c,d}}]/Det[{{a,b}, {c,d}}]]
Inverse[{{1,2}, {2,3}}] . {{1,2}, {2,3}}
[,1] [,2]
[1,] 1 2
[2,] 3 4
[1] TRUE
[1] FALSE
[1] 5.3722813 -0.3722813
[1] 5.4649857 0.3659662
[,1] [,2]
[1,] -2.0 1.0
[2,] 1.5 -0.5
[,1] [,2]
[1,] -2.0 1.0
[2,] 1.5 -0.5
[,1] [,2]
[1,] -2.0 1.0
[2,] 1.5 -0.5
Scilab 2023.1.0 (May 23 2023, 09:23:00)
ans =
-0.3722813 + 0.i
5.3722813 + 0.i
ans =
5.4649857
0.3659662
ans =
-2. 1.
1.5 -0.5
Se a matriz \(A\) é quadrada e singular, isto é, se \(\text{rank}(A)<k\) ou \(\det(A)=0\), então existe pelo menos matriz inversa denominada pseudoinversa ou inversa generalizada que é denotada por \(A^{-}\) tal que \(A=AA^{-}A\).
Moore–Penrose inverse
: Wikipedia
[,1] [,2]
[1,] 1 1
[2,] 2 2
[1] FALSE
Error in solve.default(x) :
Lapack routine dgesv: system is exactly singular: U[2,2] = 0
[,1] [,2]
[1,] 0.1 0.2
[2,] 0.1 0.2
[,1] [,2]
[1,] 1 1
[2,] 2 2
[1] TRUE
[,1] [,2]
[1,] 0.1 0.2
[2,] 0.1 0.2
Scilab 2023.1.0 (May 23 2023, 09:23:00)
x =
1. 1.
2. 2.
xi =
0.1 0.2
0.1 0.2
ans =
1. 1.
2. 2.
O índice da matriz quadrada \(A\) é a quantidade de autovalores positivos, isto é, \(\text{index}(A)=\#\{\lambda_i>0\}_{i=1}^{k}\), sendo \(\#\) o operador de contagem dos autovalores positivos do espectro ordenado de \(A\).
A assinatura da matriz quadrada \(A\) é a quantidade de autovalores positivos menos os negativos, isto é, \(\text{signature}(A)=\#\{\lambda_i>0\}_{i=1}^{k}-\#\{\lambda_i<0\}_{i=1}^{k}\).
Assim, a assinatura de uma matriz quadrada real é dada por \(\text{signature}(A)=2\times \text{index}(A)-\text{rank}(A)\).
A matriz quadrada real \(\underset{k\times k}{A}\) pode ser classificada como:
Se a matriz \(A\) é (semi)definida positiva, então \(-A\) é (semi)definida negativa.
PositiveDefiniteMatrixQ[{{1,0,0},{0,2,0},{0,0,3}}]
PositiveDefiniteMatrixQ[{{1,0,0},{0,2,0},{0,0,0}}]
NegativeDefiniteMatrixQ[{{-1,0,0},{0,-2,0},{0,0,-3}}]
NegativeDefiniteMatrixQ[{{-1,0,0},{0,-2,0},{0,0,0}}]
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 3
[1] TRUE
[1] TRUE
[1] TRUE
[1] 6
[1] 3 2 1
[1] 3 2 1
[1] TRUE
[1] 3
[1] 3
[1] 3
[1] 3
[1] FALSE
[1] TRUE
Positive definite
[1] TRUE
Positive semi definite
[1] FALSE
if((index==0 & rank==dim) | signature==-dim) {cat("Negative definite")}
matrixcalc::is.negative.semi.definite(x)
[1] FALSE
if(index==0 | signature==-rank) {cat("Negative semi definite")}
x <- matrix(c(-1, 0, 0,
0, -2, 0,
0, 0, -3),
nrow=3, ncol=3, byrow=TRUE)
x
[,1] [,2] [,3]
[1,] -1 0 0
[2,] 0 -2 0
[3,] 0 0 -3
[1] TRUE
[1] TRUE
[1] TRUE
[1] -6
[1] -1 -2 -3
[1] 3 2 1
[1] FALSE
[1] 0
[1] 3
[1] -3
[1] 3
[1] FALSE
[1] FALSE
if(dim==rank & rank==index & index==signature) {cat("Positive definite")}
matrixcalc::is.positive.semi.definite(x)
[1] FALSE
if(rank==index & index==signature) {cat("Positive semi definite")}
matrixcalc::is.negative.definite(x)
[1] TRUE
Negative definite
[1] TRUE
Negative semi definite
[,1] [,2] [,3]
[1,] -1 0 0
[2,] 0 -2 0
[3,] 0 0 0
[1] TRUE
[1] TRUE
[1] FALSE
[1] 0
[1] 0 -1 -2
[1] 2 1 0
[1] FALSE
[1] 0
[1] 2
[1] -2
[1] 3
[1] FALSE
[1] FALSE
if(dim==rank & rank==index & index==signature) {cat("Positive definite")}
matrixcalc::is.positive.semi.definite(x)
[1] FALSE
if(rank==index & index==signature) {cat("Positive semi definite")}
matrixcalc::is.negative.definite(x)
[1] FALSE
if((index==0 & rank==dim) | signature==-dim) {cat("Negative definite")}
matrixcalc::is.negative.semi.definite(x)
[1] TRUE
Negative semi definite
[,1] [,2] [,3]
[1,] 1 1 1.0
[2,] 1 1 1.0
[3,] 1 1 0.5
[1] TRUE
[1] TRUE
[1] FALSE
[1] 0
[1] 2.850781e+00 -9.436896e-16 -3.507811e-01
[1] 2.8507811 0.3507811 0.0000000
[1] FALSE
[1] 1
[1] 2
[1] 0
[1] 3
[1] TRUE
Indefinite
[1] FALSE
if(dim==rank & rank==index & index==signature) {cat("Positive definite")}
matrixcalc::is.positive.semi.definite(x)
[1] FALSE
if(rank==index & index==signature) {cat("Positive semi definite")}
matrixcalc::is.negative.definite(x)
[1] FALSE
if((index==0 & rank==dim) | signature==-dim) {cat("Negative definite")}
matrixcalc::is.negative.semi.definite(x)
[1] FALSE
O determinante real de uma matriz simétrica não singular \(A\) é dado pelo produtório dos autovalores, isto é, \(\det(A)=\prod_{i=1}^{k}{\lambda_i}\), sendo que o valor do determinante pertence ao conjunto dos reais.
A matriz quadrada, real, simétrica e não singular \(A\) não é necessariamente definida positiva, pois seu espectro ordenado não é necessariamente positivo.
Uma matriz quadrada, real e simétrica \(A\) é definida positiva quando todos os seus autovalores são maiores que zero:
\[ \lambda_i > 0 \quad \text{para todo } i \]
Se algum autovalor for zero, dizemos que \(A\) é semidefinida positiva, não definida positiva.
Por exemplo:
\[ A = \begin{pmatrix} 2 & 0 \\ 0 & 3 \end{pmatrix} \quad \Rightarrow \quad \text{autovalores: } 2, 3 \quad \Rightarrow \text{definida positiva} \]
\[ B = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \quad \Rightarrow \quad \text{autovalores: } 1, 0 \quad \Rightarrow \text{semidefinida positiva} \]
\[ C = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \quad \Rightarrow \quad \text{autovalores: } \pm 1 \quad \Rightarrow \text{indefinida} \]
Logo, se a matriz é simétrica e todos os autovalores são maiores que zero, ela é definida positiva. Se for simétrica, mas tiver pelo menos um autovalor negativo ou nulo, não é.
A média geométrica dos autovalores é \(\sqrt[k]{\det(A)}=\sqrt[k]{\prod_{i=1}^{k}{\lambda_i}}\).
Se a matriz é definida positiva, então seu determinando é positivo.
O espectro ordenado de matriz quadrada expressa uma parte importante, mas não completa, de seu conteúdo informacional.
Uma matriz quadrada definida positiva \(A\) é denotada matematicamente por \(A>0\).
Uma matriz quadrada é ortogonal se a inversa é igual à transposta, i.e., \(A^{-1}=A^{\prime}\). Portanto, \(AA^{\prime}=I\).
O valor absoluto do determinante de matriz ortogonal é igual a 1, i.e., \(|\det(A)|=1\).
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 3
[1] TRUE
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 2 0
[3,] 0 0 3
[1] FALSE
[,1] [,2] [,3] [,4]
[1,] 0 0 0 1
[2,] 0 0 1 0
[3,] 1 0 0 0
[4,] 0 1 0 0
[1] TRUE
orthogonal matrix {{1/sqrt(2),-1/sqrt(2)},{1/sqrt(2),1/sqrt(2)}}
orthogonal matrix {{cos(x), -sin(x)}, {sin(x), cos(x)}}
Transpose[{{cos(x), -sin(x)}, {sin(x), cos(x)}}] . {{cos(x), -sin(x)}, {sin(x), cos(x)}} // Simplify
A soma dos valores da diagonal principal de uma matriz quadrada real é denominada traço da matriz \(A\) e é denotado por \(\mathrm{tr}(A)=\sum_{i=1}^{k}{a_{i,i}}\).
O traço de \(A\) também pode ser expresso pelo somatório de seus autovalores, isto é, \(\mathrm{tr}(A)=\sum_{i=1}^{k}{\lambda_i}\).
A média aritmética dos autovalores é obtida dividindo o traço por \(k\), isto é, \(\frac{\mathrm{tr}(A)}{k}=\frac{1}{k}\sum_{i=1}^{k}{\lambda_i}\). Se o espectro ordenado de \(A\) é definida positiva, então o traço é positivo.
O traço da matriz identidade de dimensão \(k\times k\) é \(k\), isto é, \(\mathrm{tr}(I_{k,k})=k\).
1/2 (a + d - sqrt(a^2 + 4 b c - 2 a d + d^2)) + 1/2 (a + d + sqrt(a^2 + 4 b c - 2 a d + d^2))
[,1] [,2]
[1,] 1 -5
[2,] -5 1
[1] 2
[1] 6 -4
Scilab 2023.1.0 (May 23 2023, 09:23:00)
x =
1. -5.
-5. 1.
ans =
2.
ans =
-4.
6.
O espectro ordenado da matriz expressa uma parte importante de seu conteúdo informacional.
No entanto, apenas os autovalores não expressam toda a informação de uma matriz. Os autovalores e seus respectivos autovetores expressam toda a informação de uma matriz quadrada e não singular.
Um método de fatoração de matriz simétrica definida positiva (SDP) é a decomposição de Cholesky.
A decomposição de Cholesky extrai uma matriz triangular incluindo a diagonal principal.
A matriz \(L\) é matriz triangular de Cholesky:
\[ \begin{align} L&=\left[\begin{array}{rr} \sqrt{2} & \frac{1}{\sqrt{2}}\\ 0 & \sqrt{\frac{3}{2}} \end{array}\right]\\ \left[\begin{array}{rr} 2 & 1\\ 1 & 2 \end{array}\right]&=L^{\prime}L \end{align} \]
Outro método de fatoração é a decomposição espectral. A matriz de autovetores ortonormalizados é ortogonal.
A matriz \(E\) é o produto da matriz de autovetores ortonormalizados pelo vetor da raiz quadrada dos autovalores:
\[ \begin{align} E&=\left[\begin{array}{rr} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{array}\right] \left[\begin{array}{r} 1\\ \sqrt{3} \end{array}\right]\\ \left[\begin{array}{rr} 2 & 1\\ 1 & 2 \end{array}\right]&=E^{\prime}E \end{align} \]
PositiveDefiniteMatrixQ{{2,1},{1,2}}]
CholeskyDecomposition[{{2,1},{1,2}}]
Transpose[{{sqrt(2), sqrt(2)/2}, {0, sqrt(3/2)}}] . {{sqrt(2), sqrt(2)/2}, {0, sqrt(3/2)}}
DiagonalMatrix[Sort[Eigenvalues[{{2,1},{1,2}}]]]
Orthogonalize[Eigenvectors[{{2,1},{1,2}}]]
[,1] [,2]
[1,] 2 1
[2,] 1 2
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[,1] [,2]
[1,] 1.414214 0.7071068
[2,] 0.000000 1.2247449
[1] TRUE
[,1] [,2]
[1,] 0.7071068 0.7071068
[2,] -0.7071068 0.7071068
[,1] [,2]
[1,] 1 0.000000
[2,] 0 1.732051
[,1] [,2]
[1,] 0.7071068 1.224745
[2,] -0.7071068 1.224745
[1] TRUE
[,1] [,2]
[1,] 1 0
[2,] 0 3
[1] TRUE
Scilab 2023.1.0 (May 23 2023, 09:23:00)
A =
2. 1.
1. 2.
L =
1.4142136 0.7071068
0. 1.2247449
ans =
2. 1.
1. 2.
v =
-0.7071068 0.7071068
0.7071068 0.7071068
e =
1. 0.
0. 3.
ans =
2. 1.
1. 2.
E =
-0.7071068 1.2247449
0.7071068 1.2247449
ans =
2. 1.
1. 2.
Alguns resultados sobre matrizes reais, quadradas e não singulares de mesma dimensão envolvendo inversão, determinante e traço são importantes para alguns cálculos:
A ideia central é identificar quando um conjunto de vetores (ou equações) não traz informação nova porque um ou mais deles podem ser expressos como combinações dos outros.
No exemplo dado:
\[ \begin{aligned} 3x - 2y &= 7 \\ -6x + 4y &= -14 \end{aligned} \]
A segunda equação é o dobro da primeira multiplicada por \(-2\). Ou seja, ela não adiciona nenhuma informação nova. Assim, dizemos que as duas equações são linearmente dependentes: uma é combinação linear da outra.
A noção de dependência linear é mais importante não apenas para equações, mas também para uma variedade de outras aplicações.
Se interpretarmos os vetores como segmentos de reta direcionados, então dois vetores dependentes têm direções iguais ou opostas, conforme mostrado na Fig. 14.19.
Fig. 14.19. Pares de vetores linearmente dependentes.
Generalizando para vetores:
Se temos vetores \(a_1, a_2, \ldots, a_n\), dizemos que eles são linearmente dependentes se existem coeficientes \(\lambda_1, \ldots, \lambda_n\), não todos nulos, tais que:
\[ \lambda_1 a_1 + \lambda_2 a_2 + \cdots + \lambda_n a_n = 0 \]
Ou, em forma matricial, \(\lambda^\top a = 0\), onde \(a = [a_1\; a_2\; \ldots\; a_n]\) e \(\lambda^\top = [\lambda_1\; \lambda_2\; \ldots\; \lambda_n]\).
Se essa equação só se satisfaz quando todos \(\lambda_i = 0\), os vetores são linearmente independentes.
Exemplo geométrico:
No caso de três vetores no espaço \(\mathbb{R}^3\):
Se
\[ \lambda_1 a_1 + \lambda_2 a_2 + \lambda_3 a_3 = 0 \quad \text{com } [\lambda_1\; \lambda_2\; \lambda_3] \ne 0, \]
então, isolando \(a_3\):
\[ a_3 = -\dfrac{\lambda_1}{\lambda_3} a_1 - \dfrac{\lambda_2}{\lambda_3} a_2 \]
Ou seja, \(a_3\) está no plano formado por \(a_1\) e \(a_2\), e não traz nova direção ao conjunto. Isso caracteriza dependência linear.
Fig. 14.20. Três vetores linearmente dependentes.
Se os três vetores forem linearmente independentes, eles apontam em direções diferentes e geram um volume tridimensional (não colapsam num plano).
A teoria geral de \(n\) equações lineares com \(n\) incógnitas contém precisamente os casos considerados nos exemplos anteriores. Escrevemos o sistema na forma matricial
\[Ax=b\]
Caso 1: \(\det(A) \ne 0\) (matriz não singular): Existe uma solução única.
Subcaso 1a: \(b = 0\) (equações homogêneas): Existe apenas a solução trivial \(x = 0\).
Subcaso 1b: \(b \ne 0\) (equações não homogêneas): Existe uma solução única que não é trivial, ou seja, pelo menos um \(x_i\) é diferente de zero.
Caso 2: \(\det(A) = 0\) (matriz singular): Os lados esquerdos da Eq. (14.8.4) são linearmente dependentes.
Subcaso 2a: As equações (lados esquerdo e direito combinados) são linearmente dependentes. Pelo menos uma equação pode ser eliminada. Se existe uma solução, então existem infinitas soluções.
Subcaso 2b: As equações (lados esquerdo e direito combinados) não são linearmente dependentes. Eles se contradizem. Não existe solução alguma.
A linguagem R é essencialmente vetorial e matricial, o que significa que suas operações básicas são definidas sobre vetores, matrizes e arrays, e não elemento a elemento como em linguagens escalarizadas.
Características vetoriais:
[1] 2 3 4
[1] TRUE
[1] 2
[1] 1
Características matriciais:
matrix
e array
têm suporte
nativo: [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
[1] TRUE
[,1] [,2] [,3]
[1,] 6 4 2
[2,] 5 3 1
[1] TRUE
[,1] [,2] [,3]
[1,] 26 16 6
[2,] 37 23 9
[3,] 48 30 12
[1] TRUE
[1] TRUE
t()
, solve()
,
eigen()
manipulam diretamente estruturas matriciais.Essa natureza permite escrever código conciso, eficiente e expressivo, ideal para estatística, álgebra linear e modelagem numérica.
[,1] [,2] [,3]
[1,] 2 1 -1
[2,] -3 -1 2
[3,] -2 1 2
[1] 8 -11 -3
\left[
\begin{array}{rrrr}
2 & 1 & -1 & 8 \\
-3 & -1 & 2 & -11 \\
-2 & 1 & 2 & -3 \\
\end{array} \right]
\left[
\begin{array}{rrrr}
2 & 1 & -1 & 8 \\
-3 & -1 & 2 & -11 \\
-2 & 1 & 2 & -3 \\
\end{array} \right]
\left[
\begin{array}{llll}
1 & 1/2 & -1/2 & 8 \\
-3/2 & -1/2 & 1 & -11 \\
-1 & 1/2 & 1 & -3 \\
\end{array} \right]
[,1] [,2] [,3] [,4]
[1,] 1/6 1/3 1/2 2/3
[2,] 5/6 1 7/6 4/3
[3,] 3/2 5/3 11/6 2
\left[
\begin{array}{rrrr}
0.16667 & 0.33333 & 0.50000 & 0.66667 \\
0.83333 & 1.00000 & 1.16667 & 1.33333 \\
1.50000 & 1.66667 & 1.83333 & 2.00000 \\
\end{array} \right]
[,1] [,2] [,3]
[1,] "a" "b" "c"
[,1] [,2] [,3]
[1,] "a" "c" "e"
[2,] "b" "d" "f"
[1] 1 3
[1] "matrix" "array"
[1] 2 3
[1] "matrix" "array"
[,1]
[1,] "a"
[2,] "b"
[3,] "c"
[1] 3 1
[,1] [,2]
[1,] "a" "b"
[2,] "c" "d"
[3,] "e" "f"
[1] 3 2
[,1]
[1,] "(a) * (a) + (c) * (b) + (e) * (c)"
[2,] "(b) * (a) + (d) * (b) + (f) * (c)"
[1] 2 1
[,1] [,2]
[1,] "(a) * (a) + (b) * (c) + (c) * (e)" "(a) * (b) + (b) * (d) + (c) * (f)"
[1] 1 2
[,1] [,2] [,3]
[1,] "(a) * (a)" "(b) * (b)" "(c) * (c)"
[,1] [,2] [,3]
[1,] "(a) * (a)" "(c) * (c)" "(e) * (e)"
[2,] "(b) * (b)" "(d) * (d)" "(f) * (f)"
[,1] [,2] [,3]
[1,] "a + a" "b + b" "c + c"
[,1] [,2] [,3]
[1,] "a + a" "c + c" "e + e"
[2,] "b + b" "d + d" "f + f"
[,1] [,2] [,3]
[1,] "1" "1" "1"
[,1] [,2] [,3]
[1,] "1" "1" "1"
[2,] "1" "1" "1"
[1] "(a) * (a) + (b) * (b) + (c) * (c)"
[1] "(a) * (a) + (b) * (b) + (c) * (c) + (d) * (d) + (e) * (e) + (f) * (f)"
[1] "(a) * (a) + (b) * (b) + (c) * (c)"
[1] "(a) * (a) + (b) * (b) + (c) * (c) + (d) * (d) + (e) * (e) + (f) * (f)"
[,1]
[1,] "(a) * (a) + (b) * (b) + (c) * (c)"
[,1] [,2]
[1,] "(a) * (a) + (c) * (c) + (e) * (e)" "(a) * (b) + (c) * (d) + (e) * (f)"
[2,] "(b) * (a) + (d) * (c) + (f) * (e)" "(b) * (b) + (d) * (d) + (f) * (f)"
[,1]
[1,] "(a) * (a) + (c) * (b) + (e) * (c)"
[2,] "(b) * (a) + (d) * (b) + (f) * (c)"
, , 1, 1
[,1] [,2] [,3]
[1,] "(a) * (a)" "(b) * (a)" "(c) * (a)"
, , 1, 2
[,1] [,2] [,3]
[1,] "(a) * (b)" "(b) * (b)" "(c) * (b)"
, , 1, 3
[,1] [,2] [,3]
[1,] "(a) * (c)" "(b) * (c)" "(c) * (c)"
, , 1, 1
[,1] [,2] [,3]
[1,] "(a) * (a)" "(c) * (a)" "(e) * (a)"
[2,] "(b) * (a)" "(d) * (a)" "(f) * (a)"
, , 2, 1
[,1] [,2] [,3]
[1,] "(a) * (b)" "(c) * (b)" "(e) * (b)"
[2,] "(b) * (b)" "(d) * (b)" "(f) * (b)"
, , 1, 2
[,1] [,2] [,3]
[1,] "(a) * (c)" "(c) * (c)" "(e) * (c)"
[2,] "(b) * (c)" "(d) * (c)" "(f) * (c)"
, , 2, 2
[,1] [,2] [,3]
[1,] "(a) * (d)" "(c) * (d)" "(e) * (d)"
[2,] "(b) * (d)" "(d) * (d)" "(f) * (d)"
, , 1, 3
[,1] [,2] [,3]
[1,] "(a) * (e)" "(c) * (e)" "(e) * (e)"
[2,] "(b) * (e)" "(d) * (e)" "(f) * (e)"
, , 2, 3
[,1] [,2] [,3]
[1,] "(a) * (f)" "(c) * (f)" "(e) * (f)"
[2,] "(b) * (f)" "(d) * (f)" "(f) * (f)"
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] "(a) * (a)" "(a) * (b)" "(a) * (c)" "(b) * (a)" "(b) * (b)" "(b) * (c)"
[,7] [,8] [,9]
[1,] "(c) * (a)" "(c) * (b)" "(c) * (c)"
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] "(a) * (a)" "(a) * (c)" "(a) * (e)" "(c) * (a)" "(c) * (c)" "(c) * (e)"
[2,] "(a) * (b)" "(a) * (d)" "(a) * (f)" "(c) * (b)" "(c) * (d)" "(c) * (f)"
[3,] "(b) * (a)" "(b) * (c)" "(b) * (e)" "(d) * (a)" "(d) * (c)" "(d) * (e)"
[4,] "(b) * (b)" "(b) * (d)" "(b) * (f)" "(d) * (b)" "(d) * (d)" "(d) * (f)"
[,7] [,8] [,9]
[1,] "(e) * (a)" "(e) * (c)" "(e) * (e)"
[2,] "(e) * (b)" "(e) * (d)" "(e) * (f)"
[3,] "(f) * (a)" "(f) * (c)" "(f) * (e)"
[4,] "(f) * (b)" "(f) * (d)" "(f) * (f)"
[,1] [,2] [,3]
[1,] "a" "c" "e"
[2,] "b" "d" "f"
[,1]
[1,] "a"
[2,] "b"
[3,] "c"
[4,] "d"
[5,] "e"
[6,] "f"
[,1] [,2]
[1,] "a" "b"
[2,] "c" "d"
[1] 2 2
[1] "matrix" "array"
[1] "(a)*((d)) + -(c)*((b))"
[,1] [,2]
[1,] "((d)) / ((a)*((d)) + -(c)*((b)))" "-((b)) / ((a)*((d)) + -(c)*((b)))"
[2,] "-((c)) / ((a)*((d)) + -(c)*((b)))" "((a)) / ((a)*((d)) + -(c)*((b)))"
[1] "a + d"
[,1] [,2]
[1,] "a" "b"
[2,] "c" "d"
[,1]
[1,] "a"
[2,] "c"
[3,] "b"
[4,] "d"
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10
[3,] 11 12 13 14 15
[4,] 16 17 18 19 20
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 2 3 4 5
# [2,] 6 7 8 9 10
# [3,] 11 12 13 14 15
# [4,] 16 17 18 19 20
y <- array(data=c(1:3,3:1),
dim=c(3,2),
dimname=list(c("r1","r2","r3"),
c("c1","c2"))) # i is a 3 by 2 index array
y
c1 c2
r1 1 3
r2 2 2
r3 3 1
# c1 c2
# r1 1 3
# r2 2 2
# r3 3 1
x <- matrix(c(8, 4, 5, 6),
ncol=2,
byrow=TRUE,
dimname=list(c("r1","r2"),
c("c1","c2")))
x
c1 c2
r1 8 4
r2 5 6
num [1:2, 1:2] 8 5 4 6
- attr(*, "dimnames")=List of 2
..$ : chr [1:2] "r1" "r2"
..$ : chr [1:2] "c1" "c2"
[1] TRUE
[1] FALSE
[1] TRUE
[1] FALSE
[1] FALSE
Error in matrixcalc::is.positive.definite(x) :
argument x is not a symmetric matrix
[1] 2 2
[1] 2
[1] 2
r1 r2
12 11
c1 c2
13 10
r1 r2
6.0 5.5
c1 c2
6.5 5.0
c1 c2
r1 16 8
r2 10 12
c1 c2
r1 64 16
r2 25 36
c1 c2
r1 64 16
r2 25 36
[1] 23
[1] 960
[1] 8 6
[1] 14
[1] 28
r1
28
[1] 2
[1] 11.582576 2.417424
[1] 11.627607 2.408062
[1] 11.627607 2.408062
r1 r2
c1 0.21 -0.14
c2 -0.18 0.29
[1,] 0.21 -0.14
[2,] -0.18 0.29
r1 r2
r1 1 0
r2 0 1
r1 r2
c1 8 5
c2 4 6
c1 c2
r1 8 4
r2 5 6
[,1]
[1,] 8
[2,] 5
[3,] 4
[4,] 6
x <- matrix(c(8, 4, 5, 6, 9, 7),
ncol=3,
byrow=TRUE,
dimname=list(c("r1","r2"),
c("c1","c2", "c3")))
x
c1 c2 c3
r1 8 4 5
r2 6 9 7
num [1:2, 1:3] 8 6 4 9 5 7
- attr(*, "dimnames")=List of 2
..$ : chr [1:2] "r1" "r2"
..$ : chr [1:3] "c1" "c2" "c3"
[1] TRUE
[1] FALSE
[1] FALSE
[1] 2 3
[1] 2
[1] 3
r1 r2
17 22
c1 c2 c3
14 13 12
r1 r2
5.666667 7.333333
c1 c2 c3
7.0 6.5 6.0
c1 c2 c3
r1 16 8 10
r2 12 18 14
c1 c2 c3
r1 64 16 25
r2 36 81 49
c1 c2 c3
r1 64 16 25
r2 36 81 49
[1] 39
[1] 60480
[1] 8 9
[1] 17
Error in determinant.matrix(x, logarithm, ...) :
'x' deve ser uma matriz quadrada
[1] 2
Error in eigen(x) : matriz não quadrada em 'eigen'
[1] 16.073159 3.557183
Error in solve.default(x) : 'a' (2 x 3) deve ser quadrada
[,1] [,2]
[1,] 0.19 -0.10
[2,] -0.12 0.14
[3,] 0.00 0.04
c1 c2 c3
r1 8 4 5
r2 6 9 7
c1 c2 c3
r1 8 4 5
r2 6 9 7
[,1]
[1,] 8
[2,] 6
[3,] 4
[4,] 9
[5,] 5
[6,] 7
x <- matrix(c(11.834, 6.947, 6.819,
6.947, 9.364, 5.091,
6.819, 5.091, 12.532),
ncol=3,
byrow=TRUE,
dimname=list(c("x1","x2","x3"),
c("x1","x2","x3")))
x
x1 x2 x3
x1 11.834 6.947 6.819
x2 6.947 9.364 5.091
x3 6.819 5.091 12.532
num [1:3, 1:3] 11.83 6.95 6.82 6.95 9.36 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:3] "x1" "x2" "x3"
..$ : chr [1:3] "x1" "x2" "x3"
[1] TRUE
[1] TRUE
[1] TRUE
[1] FALSE
[1] FALSE
[1] TRUE
[1] 3 3
[1] 3
[1] 3
x1 x2 x3
25.600 21.402 24.442
x1 x2 x3
25.600 21.402 24.442
x1 x2 x3
8.533333 7.134000 8.147333
x1 x2 x3
8.533333 7.134000 8.147333
x1 x2 x3
x1 23.668 13.894 13.638
x2 13.894 18.728 10.182
x3 13.638 10.182 25.064
x1 x2 x3
x1 140.04356 48.26081 46.49876
x2 48.26081 87.68450 25.91828
x3 46.49876 25.91828 157.05102
x1 x2 x3
x1 140.04356 48.26081 46.49876
x2 48.26081 87.68450 25.91828
x3 46.49876 25.91828 157.05102
[1] 71.444
[1] 80771013
x1 x2 x3
11.834 9.364 12.532
[1] 33.73
[1] 524.1175
[1] 3
[1] 23.971567 6.272996 3.485437
[1] 23.971567 6.272996 3.485437
[1] TRUE
x1 x2 x3
x1 0.17 -0.10 -0.05
x2 -0.10 0.19 -0.02
x3 -0.05 -0.02 0.12
x1 x2 x3
x1 1 0 0
x2 0 1 0
x3 0 0 1
x1 x2 x3
x1 11.834 6.947 6.819
x2 6.947 9.364 5.091
x3 6.819 5.091 12.532
x <- matrix(c(11.834, 6.947, 6.819,
6.947, 9.364, 5.091,
6.819, 5.091, 12.532),
ncol=3,
byrow=TRUE,
dimname=list(c("x1","x2","x3"),
c("x1","x2","x3")))
y <- c(3.5, 7.8, 9.2)
x <- cbind(x, y)
x
x1 x2 x3 y
x1 11.834 6.947 6.819 3.5
x2 6.947 9.364 5.091 7.8
x3 6.819 5.091 12.532 9.2
x1 x2 x3 y
x1 11.834 6.947 6.819 3.50
x2 6.947 9.364 5.091 7.80
x3 6.819 5.091 12.532 9.20
z 13.500 19.400 3.450 7.45
Se \(x\) é qualquer número positivo ou negativo, o quadrado de \(x\) é sempre positivo. Portanto, nenhum número real satisfaz a equação quadrática
\[ x^2 = -1 \tag{15.1.1} \]
No entanto, ninguém gosta de um resultado que afirme “isso é impossível”. Matemáticos começaram, então, a buscar um novo tipo de número. Pode-se formalmente escrever \(x = \sqrt{-1}\), mas não é possível determinar se \(\sqrt{-1}\) é maior ou menor do que um número real dado. Durante muito tempo, acreditou-se que era uma característica essencial dos números possuir uma “magnitude” com uma ordem específica. Consequentemente, \(\sqrt{-1}\) não poderia ser chamado de número.
Por outro lado, operações algébricas envolvendo \(\sqrt{-1}\) podiam ser realizadas sem dificuldades. Essa situação levou, finalmente, a um compromisso: \(\sqrt{-1}\) foi denominado um número imaginário. A primeira letra da palavra “imaginário” foi proposta como notação:
\[ i = \sqrt{-1} \tag{15.1.2} \]
Hoje, a ideia de que os números devem necessariamente ser ordenados de acordo com seu tamanho foi abandonada. Não há nada de misterioso nos números imaginários. Eles podem ser somados, subtraídos, multiplicados e divididos. Junto com os números reais, eles formam o conjunto dos números complexos. Cada número é da forma:
\[ a + bi \]
sendo que \(a\) e \(b\) são números reais.
Os números complexos não são apenas úteis para investigações matemáticas. Eles também possuem aplicações práticas imediatas em uma variedade de problemas na física e na engenharia. Além disso, eles aparecem cada vez mais frequentemente na literatura biológica. Por exemplo, podem ser encontrados em livros de biologia como os de K. S. Cole (1968), Grodins (1963), Jacquard (1974), Lotka (1956), Milhorn (1966), Rosen (1967) e Sollberger (1965).
Na engenharia, a notação \(i\) é substituída por \(j\) para evitar confusão com a notação \(i\) ou \(I\) para a corrente elétrica.
Para fornecer aos números complexos algum tipo de “realidade”, mapeamos esses números em pontos de um plano.
Primeiramente, representamos a reta dos números reais como um eixo horizontal. A partir de agora, chamamos esse eixo de eixo real e o denotamos por \(X\) (Fig. 15.1). Esse eixo corresponde ao eixo \(x\) em um sistema de coordenadas cartesianas usual. Em seguida, desenhamos um eixo vertical passando pelo ponto que representa o número \(0\) e escolhemos a mesma unidade de medida usada no eixo horizontal. Associamos esse eixo ao eixo imaginário, que denotamos por \(Y\).
A figura abaixo representa o plano complexo, também chamado de plano de Argand:
Fig. 15.1. O plano complexo, também chamado de plano de Argand
Associamos o número imaginário \(i\) ao ponto que corresponde a \(x = 0\), \(y = 1\) em um sistema de coordenadas cartesiano usual. Portanto, o eixo vertical é chamado de eixo imaginário e é denotado por \(Y\). Isso também explica por que \(i\) é conhecido como a unidade imaginária.
Podemos formar múltiplos positivos e negativos de \(i\) e associá-los aos pontos no eixo imaginário. Um número complexo como \(3 + i\) é representado por um ponto com coordenadas \(x = 3\), \(y = 1\) em um sistema de coordenadas \(xy\) (Fig. 15.1). De forma geral, dados dois números reais \(x\) e \(y\), o número complexo
\[ z = x + iy \tag{15.2.1} \]
é associado ao ponto \((x, y)\) em um sistema de coordenadas \(xy\).
Um número complexo consiste em uma parte real e uma parte imaginária. Assim, em \(z = x + iy\), chamamos \(x\) de parte real e \(y\) de parte imaginária. Escrevemos:
\[ \text{Re}(z) = x, \quad \text{Im}(z) = y \tag{15.2.2} \]
Um número real é, na verdade, um caso especial de um número complexo, a saber,
\[ z = x + i \cdot 0 \]
De maneira análoga, um número imaginário também é um caso especial, dado por
\[ z = 0 + iy \]
Para enfatizar que a parte real é nula, às vezes chamamos \(z = iy\) de número puramente imaginário.
Existe uma correspondência biunívoca entre os pares \((x, y)\) de números reais e os números complexos \(z = x + iy\), expressa simbolicamente como
\[ (x, y) \leftrightarrow z = x + iy \tag{15.2.3} \]
Fig. 15.2. Os números complexos podem ser interpretados como vetores
Quando interpretamos \((x, y)\) como um vetor (Seção 14.4), ou seja,
\[ a = \begin{pmatrix} x \\ y \end{pmatrix} \]
também obtemos uma correspondência biunívoca:
\[ a = \begin{pmatrix} x \\ y \end{pmatrix} \longleftrightarrow z = x + iy \tag{15.2.4} \]
O valor absoluto ou módulo de \(z = x + iy\) é definido por
\[ |z| = \sqrt{x^2 + y^2} \tag{15.2.5} \]
Este é um número positivo, exceto no caso particular em que \(z = 0\), onde
\[
|z| = 0
\]
A Fórmula (15.2.5) corresponde à equação (14.4.7) para vetores. Geometricamente, \(|z|\) representa o comprimento do vetor
\[ a = \begin{pmatrix} x \\ y \end{pmatrix} \]
ou a distância do ponto \(x + iy\) até a origem \(z = 0\).
Por exemplo, se \(z = -12 + 5i\), então:
\[ |z| = \sqrt{144 + 25} = 13 \]
Também é conveniente introduzir coordenadas polares para números complexos. Na Seção 5.5, consideramos as coordenadas polares \((r, \theta)\) de um ponto \(P\) no plano \(xy\).
Chamamos \(r\) de distância do ponto \(P\) até a origem \(O\), e \(\theta\) do ângulo entre o eixo positivo de \(x\) e a reta \(OP\), medido no sentido anti-horário. Se \(P\) coincide com \(O\), então \(r = 0\) e \(\theta\) não é determinado.
Sejam \(x\) e \(y\) as coordenadas cartesianas do ponto \(P\), assumindo que ambas são medidas e plotadas na mesma unidade. Então, obtemos:
\[ x = r \cos(\theta), \quad y = r \sin(\theta) \tag{15.2.6} \]
Da mesma forma, podemos introduzir as coordenadas polares de um ponto no plano complexo (Fig. 15.3). A distância \(r\) é idêntica ao valor absoluto de \(z\), ou seja,
\[ r = |z| \tag{15.2.7} \]
Fig. 15.3. Coordenadas polares no plano complexo
O ângulo polar \(\theta\), também chamado de argumento ou ângulo de fase, varia de \(0\) a \(2\pi\) radianos (ou de \(0^\circ\) a \(360^\circ\)). Para \(z = 0\), o ângulo polar \(\theta\) não é determinado.
Para um número real \(z = x \neq 0\), o ângulo polar é \(\theta = 0\) se \(x > 0\) e \(\theta = \pi\) (ou \(180^\circ\)) se \(x < 0\). Para um número puramente imaginário \(iy\), temos:
Às vezes, é conveniente usar ângulos polares negativos. Assim, um ângulo polar de \(-\frac{\pi}{2}\) significa que o ponto está localizado no lado negativo do eixo imaginário.
Usando as fórmulas (15.2.6), podemos reescrever \(z = x + iy\) na forma:
\[ z = r \left(\cos(\theta) + i \sin(\theta)\right) \tag{15.2.8} \]
Considere o número complexo \(z = -7 +
5i\), com componentes \(x = -7\)
e \(y = 5\). O ponto correspondente
está localizado no segundo quadrante.
O valor absoluto é:
\[ r = |z| = \sqrt{49 + 25} = \sqrt{74} \approx 8.602 \]
A partir da fórmula (15.2.6), obtemos:
\[ \cos(\theta) = \dfrac{x}{r} = \dfrac{-7}{8.602} \approx -0.8137, \quad \sin(\theta) = \dfrac{y}{r} = \dfrac{5}{8.602} \approx 0.5812 \]
Consultando uma tabela de funções trigonométricas, obtemos:
\[ \theta \approx 144.4^\circ \]
[1] TRUE
[1] -7-5i
# Calcular o módulo (valor absoluto)
mod_z <- Mod(z)
# Calcular a parte real e imaginária
real_z <- Re(z)
imag_z <- Im(z)
# Calcular o argumento (ângulo polar em radianos e graus)
arg_z_rad <- Arg(z)
arg_z_deg <- arg_z_rad * (180 / pi) # Converter para graus
# Representação na forma polar
r <- mod_z
theta <- arg_z_rad
polar_form <- r * exp(1i * theta) # Deve resultar no mesmo z
# Exibir os resultados
cat("Número complexo: ", z, "\n")
Número complexo: -7+5i
Módulo: 8.602325
Parte real: -7
Parte imaginária: 5
Argumento (graus): 144.4623 °
Forma polar: -7+5i
# Plotar o número complexo no plano complexo
plot(NA, xlim = c(-10, 10), ylim = c(-10, 10), xlab = "Re", ylab = "Im",
main = "Representação do Número Complexo", asp = 1)
# Adicionar eixos
abline(h = 0, v = 0, col = "gray", lty = 2)
# Adicionar vetor representando o número complexo
arrows(0, 0, real_z, imag_z, col = "blue", lwd = 2)
# Adicionar ponto na posição do número complexo
points(real_z, imag_z, pch = 19, col = "red")
# Adicionar rótulo ao número complexo
text(real_z, imag_z, labels = paste0(" ", z), pos = 4, col = "red", cex = 1.2)
// Definir o número complexo
z = -7 + 5*%i;
// Calcular o módulo (valor absoluto)
mod_z = abs(z);
// Calcular a parte real e imaginária
real_z = real(z);
imag_z = imag(z);
// Calcular o argumento (ângulo polar em radianos e graus)
arg_z_rad = atan(imag_z / real_z); // ou simplesmente arg(z)
arg_z_deg = arg_z_rad * (180 / %pi) + 180; // Converter para graus
// Representação na forma polar
r = mod_z;
theta = arg_z_rad;
polar_form = r * exp(%i * theta); // Deve resultar no mesmo z
// Exibir os resultados
disp("Número complexo: " + string(z));
disp("Módulo: " + string(mod_z));
disp("Parte real: " + string(real_z));
disp("Parte imaginária: " + string(imag_z));
disp("Argumento (graus): " + string(arg_z_deg) + "°");
disp("Forma polar: " + string(polar_form));
quit
Scilab 2023.1.0 (May 23 2023, 09:23:00)
"Número complexo: -7+%i*5"
"Módulo: 8.6023253"
"Parte real: -7"
"Parte imaginária: 5"
"Argumento (graus): 144.46232°"
"Forma polar: 7-%i*5"
O número imaginário foi criado para resolver uma equação \(z^2 = -a\), sendo que \(a \in \mathbb{R}_{+}^{*}\). Note que a unidade imaginária não pertence ao conjunto dos reais, isto é, \(i \notin \mathbb{R}\). As raízes da equação são obtidas da seguinte maneira:
\[ \begin{align} z^2 &= -a\\ \sqrt{z^2} &= \sqrt{-a}\\ z &= \pm \sqrt{a} \sqrt{-1} \end{align} \]
A unidade imaginária é definida como \(i = \sqrt{-1}\). Dessa forma, tem-se que a solução é dada por \(n = \pm \sqrt{a} \times i\).
solve z^2 = -a
Resultado: \(z = \pm \sqrt{a} \times i\)
O número imaginário é denotado por \(bi\), sendo que \(b \in \mathbb{R}\). O conjunto dos números imaginários é dado por \(\mathbb{I} = \{z = bi \mid i = \sqrt{-1}, b \in \mathbb{R} \}\).
Os números imaginários positivos, não negativos, negativos e não positivos são denotados, respectivamente, por \(\mathbb{I}_+\), \(\mathbb{I}_+^*\), \(\mathbb{I}_-\) e \(\mathbb{I}_-^*\).
Is i a pure imaginary?
Resultado: True
Is 1 a pure imaginary?
Resultado: False
Is 0 a pure imaginary?
Resultado: True
A unidade imaginária possui as seguintes propriedades:
\[ i^{1+4k} = i, \quad i^{2+4k} = -1, \quad i^{3+4k} = -i, \quad i^{4+4k} = 1 \]
sendo que \(k \in \mathbb{N}\). Dessa forma,
\[ i^{1+4k} + i^{2+4k} + i^{3+4k} + i^{4+4k} = 0 \]
Por exemplo, \(i + i^2 + i^3 + i^4 = 0\). Como \(\sum_{j=1}^{2011} i^j = -1\), pois o resto inteiro de \(\frac{2011}{4}\) é 3, a soma
\[ \sum_{j=1}^{2008} i^j = 0 \]
e
\[ \sum_{j=1}^{2011} i^j = i^{2009} + i^{2010} + i^{2011} = i^{1+2008} + i^{2+2008} + i^{3+2008} = i -1 -i = -1 \]
A cota inferior do conjunto dos números inteiros é \(-\infty\). A cota superior do conjunto dos números inteiros é \(+\infty\).
A parte imaginária do número imaginário \(bi\) é o número real \(b\), isto é, \(\Im (bi) = b\).
O conjunto dos imaginários possui a propriedade de fechamento para as operações de adição e de subtração, isto é:
Se \(a \in \mathbb{I}\) e \(b \in \mathbb{I}\), então \(a + b \in \mathbb{R}\) e \(a - b \in \mathbb{I}\).
Um número imaginário \(bi\) é maior ou igual a outro imaginário \(ci\) se a parte imaginária do primeiro é maior ou igual à parte imaginária do segundo e vice-versa, isto é, \(bi \geq ci \iff b \geq c\).
O número complexo é uma extensão natural dos números imaginários. O número complexo é necessário para solucionar equações polinomiais de grau maior ou igual a 2.
Considere a equação polinomial de grau 2 ou equação quadrática \(z^2 + z + 2 = 0\). Note que \(z \notin \mathbb{R}\), isto é, não existe solução real ou imaginária para essa equação quadrática. Dessa forma, outro tipo de número faz-se necessário para resolvê-la. As raízes da equação são obtidas da seguinte maneira:
\[ \begin{align} z^2 + z + 2 &= 0 \\ z^2 + z &= -2\\ z^2 + z + \dfrac{1}{4} &= -2 + \dfrac{1}{4}\\ \left(z + \dfrac{1}{2} \right)^2 &= -\dfrac{7}{4}\\ z + \dfrac{1}{2} &= \pm \sqrt{-\dfrac{7}{4}}\\ z &= -\dfrac{1}{2} \pm \dfrac{\sqrt{7} i}{2} \end{align} \]
solve z^2 + z + 2 = 0
Resultado: \(z = -\dfrac{1}{2} \pm \dfrac{\sqrt{7} i}{2}\)
Dessa forma, tem-se que o conjunto dos números complexos é definido
por:
\[
\mathbb{C} = \{ z = a + bi \mid i = \sqrt{-1}, a \in \mathbb{R}, b \in
\mathbb{R} \}
\]
O número \(bi\) é chamado de número complexo puro. Note que \(bi = ib\).
Os números complexos, com a ausência do número 0, é denotado por \(\mathbb{C}^*\).
A parte real do número complexo \(z = a + bi\) é denotada por \(\Re(z) = a\) e sua parte imaginária, por \(\Im(z) = a\). Se a parte imaginária do complexo \(z = a + bi\) é nula, isto é, \(\Im(z) = 0\), então \(z = a\), isto é, \(z\) é um número real. Um número complexo não real tem a parte imaginária não nula, isto é, \(z = a + bi\), \(\Im(z) = b \neq 0\).
O número complexo \(-\dfrac{1}{2} - \dfrac{\sqrt{7} i}{2}\) é chamado de conjugado do complexo \(-\dfrac{1}{2} + \dfrac{\sqrt{7} i}{2}\). De forma geral, o complexo \(\bar{z} = a - bi\) é denominado complexo conjugado de \(z = a + bi\).
Real (1 + 2i)
Resultado: 1
Imaginary (1 + 2i)
Resultado:
2
Fig. Conjuntos numéricos: natural até complexo.
Decorre da definição do número complexo que as seguintes relações entre os conjuntos de números são válidas:
\[ \mathbb{N}^* \subset \mathbb{N} \subset \mathbb{Z} \subset \mathbb{Q} \subset \mathbb{R} \subset \mathbb{C} \quad \text{e} \quad \mathbb{I} \subset \mathbb{C} \]
O conjunto dos complexos possui a propriedade de fechamento para as operações de adição, subtração, multiplicação e divisão, isto é:
Se \(a \in \mathbb{C}\) e \(b \in \mathbb{C}\), então:
Is i a complex number?
Resultado: True
Is 1+i a complex number?
Resultado: True
Is 1 a complex number?
Resultado: True
Is 0 a complex number?
Resultado: True
Um número complexo \(z = a + bi\) é igual a outro número complexo \(z^{\prime} = c + di\) se e somente se \(a = c\) e \(b = d\). Não é possível estabelecer uma relação de ordem entre os complexos. Dessa forma, o conjunto dos números complexos não possui cotas inferior e superior.
O número complexo conjugado do número complexo \(z = a + bi\) é denotado por \(\overline{z} = a - bi\). Note que o conjugado possui a mesma parte real do complexo original, isto é, \(\Re(\overline{z}) = \Re(z)\), mas a parte imaginária do conjugado é de sinal contrário à parte imaginária do complexo original, isto é, \(\Im(\overline{z}) = -\Im(z)\).
Um número complexo e seu conjugado possuem o mesmo módulo, isto é, \(|z| = |\overline{z}|\). Os números complexos conjugados possuem as seguintes propriedades:
Conjugate 1+2i
Resultado: 1 - 2i
O complexo expresso por \(z = a + bi\) é dito ter a forma retangular.
As operações de adição, subtração, multiplicação e divisão envolvendo os complexos na forma retangular \(z_1 = a + bi\) e \(z_2 = c + di\) são naturalmente expressas, respectivamente, por:
\[ \begin{align} z_1 + z_2 &= (a + bi) + (c + di) \\ z_1 + z_2&= (a + c) + (b + d)i\\\\ z_1 - z_2 &= (a + bi) - (c + di) \\ z_1 - z_2&= (a - c) + (b - d)i\\\\ z_1 \times z_2 &= (a + bi) \times (c + di) \\ &= ac + adi + cbi + bdi^2 \\ &= ac + (ad + bc)i - bd \\ z_1 \times z_2&= (ac - bd) + (ad + bc)i\\\\ \dfrac{z_1}{z_2} &= \dfrac{a + bi}{c + di} \times \dfrac{c - di}{c - di} \\ \dfrac{z_1}{z_2}&= \dfrac{ac + bd}{c^2 + d^2} + \dfrac{bc - ad}{c^2 + d^2} i, \quad z_2 \neq 0 \end{align} \]
(a+bi)+(c+di)
Resultado: a + c + (b + d)i
(a+bi)(c+di)
Resultado: ac - bd + (ad + bc) i
(a+bi)/(c+di)
Resultado: \(\dfrac{a c + b d}{c^2 + d^2} + \dfrac{ b c - a d}{c^2 + d^2}i\)
A fórmula a seguir, denominada de expansão de Maclaurin, computa o número real \(e^x\), \(x \in \mathbb{R}\), com qualquer precisão numérica:
\[ e^x = 1 + x + \dfrac{x^2}{2} + \dfrac{x^3}{2 \times 3} + \dfrac{x^4}{2 \times 3 \times 4} + \dots = \sum_{k=0}^{\infty} \dfrac{x^k}{k!}, \quad x \in \mathbb{R}, \quad n \in \mathbb{N}. \]
O resultado de \(e^z\) possui forma de número complexo, isto é, \(e^z = a + bi\), sendo que \(a \in \mathbb{R}\) e \(b \in \mathbb{R}\). Por exemplo,
\[ e^{1 + \frac{\pi}{2} i} = e i \]
O exponencial complexo possui a seguinte propriedade:
\[ e^{\overline{z}} = \overline{e^z} \]
O número complexo \(e i\) é denominado número complexo transcendete, pois a parte imaginária é transcendente, isto é, \(\Im(z) = e \in \mathbb{T}\). Se a parte real ou imaginária é transcendente, o número complexo \(z\) é transcendente, isto é, \(z = a + bi \in \mathbb{C}\) é transcendente se \(\Re(z) = a \in \mathbb{T}\) ou \(\Im(z) = b \in \mathbb{T}\). Caso contrário, o número complexo é algébrico. Isto é:
\[ \text{Complexo} \begin{cases} \text{Algébrico} \\ \text{Transcendente} \end{cases} \]
exp(1)^(1+(pi/2)i)
Resultado: ei
Is ei transcendental?
Resultado: True
Is 1 algebraic?
Resultado: True
O cômputo de \(e^{1 + \frac{\pi}{2} i}\) usando a expansão infinita anterior necessita do cômputo de infinitas potências inteiras positivas de $1 + i $, isto é, de
\[ \left(1 + \frac{\pi}{2} i \right)^2, \quad \left(1 + \frac{\pi}{2} i \right)^3, \dots \]
Para o cômputo das potências inteiras positivas de um número complexo usa-se a expansão binomial de Newton:
\[ (a + bi)^k = \sum_{j=0}^{k} \binom{k}{j} a^{k-j} (bi)^j, \quad k \in \mathbb{N}, \quad \binom{k}{j} = \frac{k!}{j! (k - j)!} \]
Os coeficientes binomiais \(\binom{k}{j}\) fixando \(k\) e variando \(j\) de 0 até \(k\) formam cada linha do triângulo de Pascal, como mostra a figura abaixo para \(k\) entre 0 e 13.
Fig. Triângulo de Pascal (ou chinês) para k até 13.
Note que \(i^0 = 1\), \(i^1 = i = \sqrt{-1}\), \(i^2 = -1\), \(i^3 = -i\), \(i^4 = 1\), \(i^5 = i, \;\ldots\). Dessa forma, tem-se que:
\[ \begin{align} \left(1 + \dfrac{\pi}{2} i \right)^k &= \sum_{j=0}^{k} \binom{k}{j} (1)^{k-j} \left(\dfrac{\pi}{2} i \right)^j \\ &= \sum_{j=0}^{k} \binom{k}{j} \left(\dfrac{\pi}{2} \right)^j i^j \end{align} \]
\[ \begin{align} \left(1 + \dfrac{\pi}{2} i \right)^k &= 1 + k \left(\dfrac{\pi}{2} \right) i + \binom{k}{2} \left(\dfrac{\pi}{2} \right)^2 (-1) + \binom{k}{3} \left(\dfrac{\pi}{2} \right)^3 (-i) \\ &\quad + \binom{k}{4} \left(\dfrac{\pi}{2} \right)^4 (1) + \binom{k}{5} \left(\dfrac{\pi}{2} \right)^5 (i) + \dots + k \left(\dfrac{\pi}{2} \right)^{k-1} i^{k-1} + \left(\dfrac{\pi}{2} \right)^k i^k \end{align} \]
Dessa forma, as potências inteiras positivas de 2 a 4 são dadas por:
\[ \begin{align} \left(1 + \dfrac{\pi}{2} i \right)^2 &= \left(1 - \dfrac{\pi^2}{4} \right) + \pi i, \\ \left(1 + \dfrac{\pi}{2} i \right)^3 &= \left(1 - \dfrac{3\pi^2}{4} \right) + \left(\dfrac{3\pi}{2} - \dfrac{\pi^3}{8} \right) i, \\ \left(1 + \dfrac{\pi}{2} i \right)^4 &= \left(1 - \dfrac{3\pi^2}{2} + \dfrac{\pi^4}{16} \right) + \left(2\pi - \dfrac{\pi^3}{2} \right) i \end{align} \]
Expand (1+(pi/2)i)^2
Resultado: \(1 - \dfrac{\pi^2}{4} + \pi i\)
Finalmente, tem-se que:
\[ \begin{align} e^{1+\frac{\pi}{2} i} &= \sum_{k=0}^{\infty} \dfrac{1}{k!} \sum_{j=0}^{k} \binom{k}{j} \left(\dfrac{\pi}{2} \right)^j i^j \\ &= 1 + \left(1 + \dfrac{\pi}{2} i \right) + \dfrac{\left(1 - \dfrac{\pi^2}{4} \right) + \pi i}{2} \\ &\quad + \dfrac{\left(1 - \dfrac{3\pi^2}{4} \right) + \left(\dfrac{3\pi}{2} - \dfrac{\pi^3}{8} \right) i}{2 \times 3} \\ &\quad + \dfrac{\left(1 - \dfrac{3\pi^2}{2} + \dfrac{\pi^4}{16} \right) + \left(2\pi - \dfrac{\pi^3}{2} \right) i}{2 \times 3 \times 4} + \cdots \end{align} \]
O resultado correto dessa soma infinita enumerável de termos pode ser obtido pelo software Mathematica 12.3 da seguinte forma:
Limit[Sum[(Sum[Binomial[k,j] ((Pi/2)^j) (I^j),{j,0,k}])/Factorial[k],{k,0,Infinity}],k->Infinity]
Resultado: ei
Outra forma de computar \(e^{1+\frac{\pi}{2} i}\) é mostrada a seguir.
O número complexo \(z = a+bi\) possui uma maneira alternativa de expressão denominada forma polar, isto é,
\[ z = r e^{\theta i} \]
sendo que \(r \in \mathbb{R}_+^*\) é a norma ou módulo, \(\theta \in \mathbb{R}\), tal que o domínio ou ramo principal é dado por \(0 < \theta < 2\pi\), é o argumento principal. Note que \(e^\theta > 0\) e que o complexo é igual a 0 apenas se \(r = 0\).
O complexo \(z\) na forma polar corresponde a um vetor de comprimento dado pelo módulo
\[ |z|=\sqrt{a^2 + b^2} \]
e pelo ângulo dado pelo argumento \(\arg(z) = \theta\). Note que o módulo do complexo é real positivo, isto é, \(|z| \in \mathbb{R}^+\) e o argumento pode ser expresso em radiano (rad) ou grau (°). O argumento em grau é dado por
\[ \theta = \dfrac{a}{\pi} \times 180^\circ \]
sendo \(a\) expresso em radiano. O argumento e o módulo são chamados de coordenadas polares.
modulus a+bi
Resultado: \(\sqrt{a^2 +
b^2}\)
|a+bi|
Resultado: \(\sqrt{a^2 +
b^2}\)
i^2
Resultado: \(-1\)
A notação polar facilita a multiplicação e divisão nos complexos, pois se \(z_1 = r_1 e^{\theta_1 i}\) e \(z_2 = r_2 e^{\theta_2 i}\), então:
\[ \begin{align} z_1 \times z_2 &= r_1 e^{\theta_1 i} \times r_2 e^{\theta_2 i} \\ z_1 \times z_2&= r_1 r_2 e^{(\theta_1 + \theta_2)i} \end{align} \]
\[ \begin{align} \dfrac{z_1}{z_2} &= \dfrac{r_1 e^{\theta_1 i}}{r_2 e^{\theta_2 i}} \\ \dfrac{z_1}{z_2}&= \dfrac{r_1}{r_2} e^{(\theta_1 - \theta_2)i}, \quad r_2 \neq 0 \end{align} \]
Por exemplo,
\[ i^{-1} = \dfrac{1}{i} = \dfrac{1}{e^{\frac{\pi}{2} i}} = e^{-\frac{\pi}{2} i} \]
Pelo fato do número complexo \(z\) ser um vetor, ele também pode ser expresso pela forma trigonométrica
\[ z = r (\cos (\theta) + \sin (\theta) i) \]
Note que as partes real e imaginária do complexo são dadas, respectivamente, por \(\Re(z) = r \cos (\theta)\) e \(\Im(z) = r \sin (\theta)\). Dessa forma, o módulo do complexo é dado por
\[ \begin{align} |z| &= \sqrt{(r \cos (\theta))^2 + (r \sin (\theta))^2} \\ &= r \sqrt{(\cos (\theta))^2 + (\sin (\theta))^2} \\ |z|&= r \end{align} \]
A seguinte identidade é, por consequência, válida:
\[ e^{\theta i} = \cos (\theta) + \sin (\theta) i \]
Como o argumento possui período \(2\pi\), tem-se que
\[ \cos (\theta) = \cos (\theta + 2\pi k) \quad \text{e} \quad \sin (\theta) = \sin (\theta + 2\pi k) \]
sendo que \(k \in \mathbb{Z}\). Dessa forma,
\[ e^{\theta i} = e^{(\theta + 2\pi k)i} \]
Portanto,
\[ z = r e^{\theta i} = r e^{(\theta + 2\pi k)i} \]
Note que a forma polar restrita ao domínio principal evita redundâncias.Fig. Valores notáveis da circunferência trigonométrica.
A fórmula de Euler é dada por
\[ z^n = r^n e^{n\theta i} \]
sendo que \(n \in \mathbb{Q}\). Ela permite computar, por exemplo, \(z^2\), \(\sqrt{z}\) e \(1/z\). Por exemplo, para \(z = i = 0 + 1i\), tem-se que
\[ z^2 = i^2 = (\sqrt{-1})^2 = -1 \]
usando apenas a definição de unidade imaginária. O módulo de \(z = i\) é 1 e o argumento \(\theta\) é \(\pi/2\). Em notação polar, tem-se
\[ i^2 = 1^2e^{2\frac{\pi}{2} i} = e^{\pi i} = \cos (\pi) + \sin (\pi) i = -1 + 0 \times i = -1 \]
|i|
Resultado: 1
argument i
Resultado: \(\pi/2\)
i^(-1) = -i
Resultado: True
|-i|
Resultado: 1
argument -i
Resultado: \(-\pi/2\)
A fórmula de De Moivre é dada por
\[ z^n = r^n (\cos (n\theta) + \sin (n\theta) i), \]
sendo que \(n \in \mathbb{Q}\).
A seguinte identidade é por consequência válida:
\[ e^{n\theta i} = \cos (n\theta) + \sin (n\theta) i \]
Portanto, por exemplo:
\[ \begin{align} z^2 &= r e^{2\theta i} \\ (0 + 1i)^2 &= \sqrt{0^2 + 1^2} \left( \cos \left( 2\dfrac{\pi}{2} \right) + \sin \left( 2\dfrac{\pi}{2} \right) i \right) \\ i^2 &= \cos (\pi) + \sin (\pi) i \\ i^2 &= -1 \end{align} \]
Note que \(e^{\pi i} = \cos (\pi) + \sin (\pi) i = -1\).
ComplexExpand[Exp[a+bI]]
Resultado: \(e^a (\cos (b) + \sin
(b)i)\)
ComplexExpand[Exp[Pi I]]
Resultado: \(-1\)
Pela fórmula de De Moivre pode-se calcular a potência racional do número complexo \(e^z\) pode ser calculada por meio da fórmula
\[ (e^z)^n = e^{zn} = e^{na} (\cos (n\theta) + \sin (n\theta) i) \]
sendo que \(z = a + \theta i\), \(a \in \mathbb{R}\) e \(n \in \mathbb{Q}\).
Sua dedução é mostrada a seguir:
\[ \begin{align} e^{\theta i} &= \cos (\theta) + \sin (\theta) i \\\\ e^{z} &= e^{a+\theta i} =e^a e^{\theta i} = e^a (\cos (\theta) + \sin (\theta) i) \\ e^{nz} &=e^{na+n\theta i} = e^{na} e^{n\theta i} = e^{na} (\cos (n\theta) + \sin (n\theta) i). \end{align} \]
Esta última fórmula permite computar, por exemplo,
\[ e^{\frac{\pi }{2}i} \]
O módulo de \(e^{\frac{\pi }{2}i}\) é 1 e o seu argumento é \(\pi/2\). Aplicando a fórmula, tem-se que:
\[ \begin{align} e^{\frac{\pi }{2}i} &= e^0 \cdot e^{\frac{\pi }{2}i} \\ &= e^0 \left( \cos \left(\dfrac{\pi}{2} \right) + \sin \left(\dfrac{\pi}{2} \right) i \right) \\ e^{\frac{\pi }{2}i}&= i \end{align} \]
Anteriormente, foi calculado o valor de \(e^{1+\frac{\pi }{2}i}\) usando a expansão de Maclaurin e o valor obtido foi \(e^{1+\frac{\pi }{2}i} = e i\). Usando a fórmula
\[ e^z = e^a (\cos (\theta) + \sin (\theta) i) \]
para computar \(e^{1+\frac{\pi i}{2}}\), tem-se que:
\[ \begin{align} e^{1+\frac{\pi }{2}i} &= e^1 \left( \cos \left(\dfrac{\pi}{2} \right) + \sin \left(\dfrac{\pi}{2} \right) i \right) \\ e^{1+\frac{\pi }{2}i}&= e i \end{align} \]
ComplexExpand[Exp[1+(Pi/2)I]]
Resultado: \(ei \approx 2.71828\; i\)
O logaritmo complexo de número complexo \(z = r e^{\theta i}\) é computado por meio da fórmula
\[ \ln (z) = \ln (r) + \theta i \]
Dessa forma, o logaritmo de número real negativo é computável.
Por exemplo,
\[ \ln (-1) = \ln (1) + \pi i = \pi i \]
ComplexExpand[Log[a+bI]]
Resultado: \(\dfrac{1}{2} \log \left(a^2 + b^2\right) + \arg (a + bi)i\)
A potência complexa do número complexo \(z = r e^{\theta i}\) pode ser calculada por meio da fórmula
\[ z^c = e^{c \ln (z)} = e^{c (\ln (r) + \theta i)} \]
sendo que \(z \in \mathbb{C}^*\) e \(c \in \mathbb{C}\). Esta fórmula permite computar, por exemplo, \(i^i\).
O módulo de \(i\) é \(r = 1\), o seu argumento é \(\theta = \pi/2\) e \(c = i\). Aplicando a fórmula, tem-se que:
\[ \begin{align} i^i &= e^{i \left(\ln (1) + \frac{\pi }{2}i\right)} \\ &= e^{\frac{\pi }{2}i^2} \\ i^i&= e^{-\frac{\pi}{2}}. \end{align} \]
Com a aplicação da fórmula da potência complexa do número complexo, obtêm-se os seguintes resultados:
\[ i^1 = 1, \quad i^0 = 1, \quad i^{2i} = e^{2i \ln i} = e^{(-i) \frac{\pi}{2}} = e^{\frac{\pi}{2}} \]
ComplexExpand[Cos[a+bI]]
Resultado: \(\dfrac{e^{-(b- ai)}+e^{(b - ai)}}{2}\)
O número complexo pode ser representado graficamente de três maneiras distintas.
A primeira representação gráfica do número complexo ocorre no plano polar \(\mathbb{R} \times \mathbb{I}\) denominado plano de Argand-Gauss, sendo que a abscissa é a reta real \(\mathbb{R}\) composta pela parte real e a ordenada é a reta imaginária \(\mathbb{I}\) composta pelo número imaginário. Por exemplo, o número complexo \(1+i\) no plano de Argand-Gauss é um vetor com origem no ponto \((0,0)\) e afixo \((1,i)\).
A segunda representação gráfica do número complexo \(z\) ocorre no plano cartesiano \(\Re(z) \times \Im(z)\) denominado plano das partes real e imaginária, sendo que a abscissa é a reta real \(\mathbb{R}\) composta pela parte real e a ordenada é a reta real \(\mathbb{R}\) composta pela parte imaginária. Por exemplo, o número complexo \(1+i\) no plano das partes real e imaginária é o afixo \((1,1)\).
A terceira representação gráfica do número complexo \(z\) ocorre no plano cartesiano \(|z| \times \arg (z)\) denominado coordenadas polares, sendo que a abscissa é a reta real \(\mathbb{R}\) composta pelo módulo e a ordenada é a reta real \(\mathbb{R}\) composta pelo argumento em radiano. Por exemplo, o número complexo \(1+i\) no plano das partes real e imaginária é o afixo
\[ \left(\sqrt{2}, \frac{\pi}{2} \right) \]
A função polinomial de grau \(n\) é definida como
\[ \begin{align} f(x)&=a_{0}+a_{1}x+a_{2}x^{2}+\cdots+a_{n-1}x^{n-1}+a_{n}x^{n} \\ f(x)&=\sum_{i=0}^{n}{a_{i}x^{i}} \end{align} \]
\(x\in\mathbb{R}\), \(a_{n}\in\mathbb{R^{\ast}}\) e \(a_{i}\in\mathbb{R}\), \(i\in \{1,2,...,n-1\}\) e \(n\in\mathbb{N}^{\ast}\).
A função polinomial de grau \(n\) possui \(D(f)=\mathbb{R}\) e \(Im(f)=\mathbb{R}\).
Se a função polinomial possui grau par e \(D(f)=\mathbb{R}\), então \(Im(f)\subset \mathbb{R}\).
Todo número real é um número complexo. Um número puramente real tem a parte imaginária nula.
Um polinômio de ordem \(n\) com coeficientes complexos possui exatamente \(n\) raízes complexas, contando com multiplicidades, de acordo com o Teorema Fundamental da Álgebra.
Se os coeficientes são complexos, as raízes não precisam ser conjugadas, diferentemente do caso em que os coeficientes são reais (onde as raízes complexas ocorrem aos pares conjugados).
Um polinômio de ordem \(n\) com coeficientes puramente reais tem \(n\) raízes complexas. As raízes complexas sempre aparecem aos pares conjugados. E.g.: Um polinômio de grau 4 com coeficientes reais pode ter: (i) 4 raízes reais; (ii) 2 raízes reais e 2 raízes complexas conjugadas; (iii) 0 raízes reais e 4 raízes complexas conjugadas (dois pares conjugados).
[1] "5*x^4 + 10*x^3 - 1*x^2 + 5*x - 1"
# A tibble: 4 × 2
root mult
<cpl> <dbl>
1 0.19187896+0.0000000i 1
2 0.05096963-0.6721664i 1
3 0.05096963+0.6721664i 1
4 -2.29381822+0.0000000i 1
[1] -2.29381822+0.0000000i 0.05096963+0.6721664i 0.05096963-0.6721664i
[4] 0.19187896+0.0000000i
[1] 0+0i 0+0i 0+0i 0+0i
[1] -2.2938144 0.1918704
[1] 165
curve(f, from = -2.5, to = 1.5, col = "blue", lwd = 2,
xlab = "x", ylab = "y",
main = expression(y == 5*x^4 + 10*x^3 - x^2 + 5*x - 1))
abline(h=0, v=0, lty=2)
# Coeficientes do polinômio em ordem decrescente de potência
coef <- c(5, 10, -1, 5, -1)
polin_strg <- pracma::poly2str(coef)
print(polin_strg)
[1] "5*x^4 + 10*x^3 - 1*x^2 + 5*x - 1"
root mult
1 0.192+0.000i 1
2 0.051-0.672i 1
3 0.051+0.672i 1
4 -2.294+0.000i 1
[1] -0.436+2.07e-25i 0.112+1.48e+00i 0.112-1.48e+00i 5.212+0.00e+00i
[1] -2.29381822+0.0000000i 0.05096963+0.6721664i 0.05096963-0.6721664i
[4] 0.19187896+0.0000000i
[1] 0+0i 0+0i 0+0i 0+0i
[1] -2.2938144 0.1918704
Compute \(\sqrt{-15 - 8i}\).
Resposta:
\[ \sqrt{-15 - 8i} = \sqrt{-1} \sqrt{15 + 8i} = i \sqrt{15 + 8i} = i (4 + i) = -1 + 4i \]
Compute
Resposta:
Usar a expansão binomial de Newton; usar o comando
Binomial[20,j]
do WolframAlpha para calcular os
coeficientes binomiais.
\[ \sum_{j=0}^{20} \binom{20}{j} 0,01^j = 1,220188 \dots \]
\[ 2^{20} = (1+1)^{20} = \sum_{j=0}^{20} \binom{20}{j} j \]
\[ (1+i)^{20} = \sum_{k=0}^{20} \binom{20}{k} i^k = -1024 \]
\[ (1 - i)^{20} = \sum_{k=0}^{20} \binom{20}{k} (-i)^k = -1024 \]
Note que
\[ (1 + i)^{4 \times n} = (1 - i)^{4 \times n} \in \mathbb{R}, \quad \text{sendo que } n \in \mathbb{N} \]
Explique a falácia decorrente da sequência de identidades
A falácia está na manipulação incorreta das propriedades da raiz quadrada. Seguindo a sequência de passos:
\[ -1 = \sqrt{-1} \cdot \sqrt{-1} = \sqrt{(-1)(-1)} = \sqrt{1} = 1 \]
O erro ocorre porque a propriedade
\[ \sqrt{a} \sqrt{b} = \sqrt{a b} \]
só é válida para números reais não negativos (\(a, b \geq 0\)). Como \(-1\) é negativo, essa manipulação não é válida.
A correção correta da identidade seria:
\[ -1 = \sqrt{-1} \sqrt{-1} \neq \sqrt{(-1)(-1)} = \sqrt{1} = 1 \]
Ou seja, \(\sqrt{-1}\) não pode ser tratado simplesmente como um número real, pois ele representa um símbolo específico.
Qual o lugar geométrico no plano cartesiano formado pelas partes real e imaginária do número complexo \(z = a + bi\) que satisfaz a equação
\[ |z + k||z - k| = k^2, \quad \text{sendo } k \in \mathbb{R}^*? \]
Resposta:
O comando
ContourPlot[Abs[a + b I - 1] Abs[a + b I + 1] == 1, {a, -2, 2}, {b, -2, 2}]
resulta na leminiscata.
ContourPlot[Abs[a + b I - 1] Abs[a + b I + 1] == 1, {a, -2, 2}, {b, -2, 2}]
ContourPlot[Abs[a + b I - 2] Abs[a + b I + 2] == 4, {a, -4, 4}, {b, -2, 2}]
Mostre que a equação
\[ x^4 + x^2 + 1 = 0 \]
é satisfeita para cada um dos seguintes complexos e seus respectivos conjugados:
\[ \dfrac{1}{2} (1 - \sqrt{3} i), \quad -\dfrac{1}{2} (1 + \sqrt{3} i) \]
Resposta:
Os conjugados de
\[ \dfrac{1}{2}(1 - \sqrt{3} i) \]
e
\[ - \dfrac{1}{2} (1 + \sqrt{3} i) \]
são, respectivamente,
\[ \dfrac{1}{2} (1 + \sqrt{3} i) \quad \text{e} \quad -\dfrac{1}{2} (1 - \sqrt{3} i) \]
Verificação da equação para
\[ x = \dfrac{1}{2} (1 - \sqrt{3} i) \]
O comando abaixo resulta TRUE
.
((1/2)*(1-sqrt(3)i))^4+((1/2)(1-sqrt(3)*i))^2+1=0
WolframAlpha: x^4+x^2+1=0
resulta as 4 raízes
complexas.
Mostre que
\[ (a + bi)(c + di)(b + ai)(d + ci) \]
sendo que as partes reais e imaginárias são não nulas, é um número real negativo.
Resposta:
\[ -(a^2 + b^2) (c^2 + d^2) \]
Compute o determinante de
\[ \left[ \begin{matrix} i & 0 & 0 \\ 0 & i & 0 \\ 0 & 0 & i \end{matrix} \right]^9 \]
Resposta:
\[ - i \]
pois
\[ \left[ \begin{matrix} i & 0 & 0 \\ 0 & i & 0 \\ 0 & 0 & i \end{matrix} \right]^9 = \left[ \begin{matrix} i & 0 & 0 \\ 0 & i & 0 \\ 0 & 0 & i \end{matrix} \right] \]
e
\[ \det \left[ \begin{matrix} i & 0 & 0 \\ 0 & i & 0 \\ 0 & 0 & i \end{matrix} \right] = i^3 = -i \]
Considere \(x = a + bi\) um número complexo.
Resposta:
\[ xi = (a + bi)i = -b + ai \]
O módulo de \(x\) e \(xi\) é:
\[ |x| = |xi| = \sqrt{a^2 + b^2} \]
As partes real e imaginária são:
\[ \Re(x) = \Im(xi) = a, \quad \Im(x) = -\Re(xi) = b \]
Por exemplo, se
\[ x = 1 + i \]
então
\[ xi = -1 + i \]
Considere a equação quadrática abaixo:
\[ z^2 - 4z + 13 = 0 \]
Solução:
Aplicamos a fórmula de Bhaskara para encontrar as raízes:
\[ z = \dfrac{-(-4) \pm \sqrt{(-4)^2 - 4(1)(13)}}{2(1)} \]
Calculando o discriminante:
\[ \Delta = (-4)^2 - 4(1)(13) = 16 - 52 = -36 \]
Como o discriminante é negativo, as raízes serão complexas conjugadas:
\[ z = \dfrac{4 \pm \sqrt{-36}}{2} = \dfrac{4 \pm 6i}{2} = 2 \pm 3i \]
Portanto, as raízes são \(z = 2 + 3i\) e \(z = 2 - 3i\), que são conjugadas complexas.
Solução geral:
A equação algébrica possui como incógnita a variável \(z\), sendo que os coeficientes são conhecidos.
A equação algébrica quadrática é dada por:
\[ a_2 z^2 + a_1 z + a_0 = 0 \]
sendo que \(z \in \mathbb{C}\) e \(a_2 \in \mathbb{C}^*\), \(a_1 \in \mathbb{C}\) e \(a_0 \in \mathbb{C}\).
As duas raízes são determinadas pela fórmula de Bhaskara:
\[ \dfrac{-a_1 \pm \sqrt{a_1^2 - 4a_0 a_2}}{2a_2}, \quad a_2\neq0 \]
O radicando \(a_1^2 - 4a_2 a_0\) é denominado discriminante. Se o discriminante é maior que zero, então as duas raízes são reais e diferentes. Se o discriminante é igual a zero, então as duas raízes são reais e iguais. Se o discriminante é menor que zero, então as duas raízes são complexas.
a2 z^2 + a1 z + a0 = 0
Resultado:
\[ z = \dfrac{-a_1 \pm \sqrt{a_1^2 - 4a_0 a_2}}{2 a_2} \quad (\text{com } a_2 \neq 0) \]
Considere a equação cúbica abaixo:
\[ 2z^3+12z^2+98z-260=0 \]
Resposta:
As raízes da equação são:
\[ z = 2, \quad z = -4 + 7i, \quad x = -4 - 7i \]
Solução geral:
A equação algébrica cúbica é dada por:
\[ a_3 z^3 + a_2 z^2 + a_1 z + a_0 = 0 \]
sendo que \(z \in \mathbb{C}\) e \(a_3 \in \mathbb{C}^*\), \(a_j \in \mathbb{C}\) e \(j \in \{0,1,2\}\).
O primeiro passo consiste em dividir os dois lados da equação algébrica cúbica por \(a_3\) e, depois, fazer a substituição da variável \(z\) por
\[ x = \dfrac{z - a_2}{3 a_3} \]
resultando assim na seguinte equação algébrica cúbica reduzida em \(x\):
\[ x^3 + r x^2 + s = 0 \]
Sendo que \(x \in \mathbb{C}\),
\[ r = \dfrac{3 a_3 a_1 - a_2^2}{3 a_3} \quad \text{e} \quad s = \dfrac{2 a_2}{27 a_3} - \dfrac{a_2 a_1}{3 a_3^2} + \dfrac{a_0}{3 a_3} \]
O discriminante é
\[ \left( \dfrac{r}{3} \right)^3 + \left( \dfrac{s}{2} \right)^2 \]
Se o discriminante é maior que zero, então uma raiz é real e duas são complexas conjugadas.
\[ x_1 = u + v \]
\[ x_{2,3} = -\dfrac{u + v}{2} \pm \dfrac{\sqrt{3}}{2} (u - v) i \]
Sendo que
\[ u = \sqrt[3]{-\dfrac{s}{2} + \sqrt{\left( \dfrac{r}{3} \right)^3 + \left( \dfrac{s}{2} \right)^2}} \quad \text{e} \quad v = \sqrt[3]{-\dfrac{s}{2} - \sqrt{\left( \dfrac{r}{3} \right)^3 + \left( \dfrac{s}{2} \right)^2}} \]
Se o discriminante é maior ou igual a zero, então as três raízes são reais e uma delas tem multiplicidade 2:
\[ x_1 = 2 \sqrt[3]{-\dfrac{s}{2}} \]
\[ x_{2,3} = -3 \sqrt[3]{-\dfrac{s}{2}} \]
Se o discriminante é menor que zero, então as três raízes são reais e diferentes:
\[ x_{1,2,3} = 2 \sqrt{-\dfrac{r}{3}} \left( \dfrac{1}{3} \arccos \left( -\dfrac{s}{2} \bigg/ \sqrt{-\dfrac{r}{3}} \right) + \dfrac{4}{3\pi} k \right) \]
sendo que \(k \in \{1,2,3\}\).
Destarte, as raízes da equação algébrica cúbica são dadas por:
\[ z_1 = x_1 - \dfrac{a_2}{3a_3}, \quad z_2 = x_2 - \dfrac{a_2}{3a_3}, \quad z_3 = x_3 - \dfrac{a_2}{3a_3} \]
Por exemplo, a equação algébrica cúbica
\[ 2z^3 + 12z^2 + 98z - 260 = 0 \]
tem a forma reduzida
\[ x^3 + 37x - 212 = 0 \]
por meio da divisão por 2 e, em seguida, pela substituição de \(z\) por
\[ x = \dfrac{z - 12}{(3 \times 4)} = x - 2 \]
As raízes da forma reduzida são
\[ 4, \quad -2 \pm 7i \]
Dessa maneira, as raízes são
\[ 4 - 2 = 2, \quad -2 \pm 7i - 2 = -4 \pm 7i \]
2z^3 + 12z^2 + 98z - 260 = 0
Resultado: \(x = -4 - 7i\) e \(x = -4 + 7i\) e \(x = 2\).
Collect (2z^3 + 12z^2 + 98z - 260)/2
Resultado: \(z^3 + 6z^2 + 49z - 130\).
Collect (x - 2)^3 + 6(x - 2)^2 + 49(x - 2) - 130
Resultado: \(x^3 + 37x - 212\).
x^3 + 37x - 212 = 0
Resultado: \(x = -2 - 7i\) e \(x = -2 + 7i\) e \(x = 4\).
As equações algébricas da forma
\[ a_2 z^{2n} + a_1 z^n + a_0 = 0 \]
e
\[ a_3 z^{3n} + a_2 z^{2n} + a_1 z^n + a_0 = 0 \]
são redutíveis respectivamente às equações algébricas quadrática e cúbica por meio das mudanças de variáveis respectivas de \(z^n\) por \(x\). As raízes são da forma
\[ z = \pm \sqrt[n]{x} \]
Por exemplo, a equação algébrica de grau quatro
\[ z^4 + 6z^2 - 16 = 0 \]
reduz-se a
\[ x^2 + 6x - 16 = 0 \]
e as duas raízes da reduzida são \(-8\) e \(2\), então as quatro raízes da equação original são
\[ \pm \sqrt{-8} = \pm 2\sqrt{2}i \quad \text{e} \quad \pm \sqrt{2} \]
A maioria, talvez até todos, os vegetais, animais, microrganismos, órgãos individuais ou partes de tecidos participam de processos cíclicos conhecidos como ritmos biológicos. Entre os ritmos biológicos mais conhecidos estão as ondas cerebrais, batimentos cardíacos, respiração, variações diárias na atividade dos rins e do fígado, movimentos diários das folhas, o ciclo mensal da menstruação e os padrões anuais de reprodução. Para uma visão geral sobre ritmos biológicos e referências, consulte Sollberger (1965).
Os processos cíclicos podem ser iniciados por sinais e forças externas, como luz, variações de temperatura e umidade, estímulos químicos e radiação corpuscular. Muitas vezes, os ritmos biológicos continuam mesmo quando os sinais ou forças externas são removidos. Por exemplo, folhas continuam seu movimento de subida e descida mesmo quando mantidas em escuridão constante. Da mesma forma, hamsters não interrompem seu ritmo diário de atividade e descanso quando mantidos sob luz constante. Isso sugere que os organismos possuem um sistema interno de temporização, conhecido como relógio biológico.
Pouco se sabe sobre a natureza dos relógios biológicos. Uma teoria sugere que esse relógio leria e contaria sinais rítmicos do ambiente físico e utilizaria essas informações para controlar o organismo. Em princípio, esse relógio funcionaria como um relógio elétrico, movido pelos sessenta ciclos por segundo da corrente alternada (cf. F. A. Brown, 1965). Outra teoria afirma que o relógio biológico é mais autônomo. Em princípio, ele se assemelharia a um relógio comum que contém seu próprio sistema oscilatório, permitindo seu funcionamento sem a necessidade de sinais periódicos externos (ver, por exemplo, Pittendrigh e Bruce, 1957).
Nesta seção, concentraremos nossa atenção na matemática dos sistemas oscilatórios. Existem diversos tipos de sistemas oscilatórios, que chamaremos de osciladores a partir de agora.
Na mecânica, conhecemos o pêndulo, o balanço de um relógio de pulso, cordas e membranas vibrantes, o diapasão e outros osciladores. No campo da eletricidade, osciladores incluem transmissores e receptores de radar e rádio.
As oscilações podem ocorrer em espaços muito pequenos, como em diversos osciladores moleculares e atômicos. Menos conhecidos são os osciladores térmicos, cujo exemplo familiar é um gêiser. Além disso, existem vários sistemas químicos capazes de gerar oscilações (para referências, ver Sollberger, 1965, p. 59).
Como exemplo instrutivo, estudaremos as oscilações de um corpo preso a uma mola (Fig. 15.9). Permitimos que o corpo se mova apenas na direção vertical. Consideramos a linha vertical como um eixo \(x\), com o lado positivo apontando para baixo e a origem \(x = 0\) localizada no ponto de equilíbrio do corpo, onde o peso do corpo e a força elástica da mola se cancelam.
Seja \(F = F(x)\) a força resultante como uma função de \(x\). Se o corpo estiver abaixo do ponto de equilíbrio, digamos em \(x_1\), a força \(F(x_1) = F_1\) será negativa, pois acelera o corpo na direção negativa de \(x\). Por outro lado, se o corpo estiver acima do ponto de equilíbrio, digamos em \(x_2\), a força \(F(x_2) = F_2\) será positiva. Portanto, \(x\) e \(F(x)\) possuem sinais opostos.
Para deslocamentos moderados de \(x\) em relação ao equilíbrio, é um fato experimental (lei de Hooke) que \(F(x)\) é uma função linear de \(x\). Assim, obtemos
\[ F(x) = -k x \tag{15.6.1} \]
sendo que \(k\) denota uma constante positiva.
Fig. 15.9. Um corpo suspenso em uma mola helicoidal como exemplo de um oscilador mecânico. A força \(F\) atuando sobre o corpo é a resultante do peso do corpo e da força elástica da mola.
O deslocamento \(x\) é uma função do tempo. Assim, podemos escrever \(x = x(t)\). Agora consideramos a velocidade \(v = v(t) = \dfrac{dx}{dt}\) e a aceleração \(\dfrac{dv}{dt} = \dfrac{d^2 x}{dt^2}\) do corpo. A segunda lei da dinâmica de Newton estabelece que
\[ M \dfrac{d^2 x}{dt^2} = F(x) \tag{15.6.2} \]
sendo que \(M\) denota a massa do corpo.
A partir de (15.6.2) e (15.6.1), segue que
\[ \boxed{M \dfrac{d^2 x}{dt^2} + kx = 0}\quad (M > 0, k > 0) \tag{15.6.3} \]
Temos assumido tacitamente que não há outras forças presentes. Em particular, excluímos o atrito e consideramos que a massa da mola é negligenciável.
A Fórmula (15.6.3) é uma equação diferencial de segunda ordem para a função desconhecida \(x = x(t)\). A equação é linear, e os coeficientes são constantes.
Esse tipo de equação diferencial é resolvido pelo método das funções de teste, que já aplicamos na Seção 11.7. Tentamos a função exponencial
\[ x = A e^{\lambda t} \tag{15.6.4} \]
Substituindo \(x\) na equação (15.6.3) por essa função, obtemos
\[ M \lambda^2 A e^{\lambda t} + k A e^{\lambda t} = (M \lambda^2 + k) A e^{\lambda t} = 0 \]
A nova equação é válida se, e somente se, \(A = 0\) ou
\[ \boxed{M \lambda^2 + k = 0} \tag{15.6.5} \]
Desconsideramos o caso trivial \(A = 0\). A equação quadrática (15.6.5) é chamada de equação característica. Ela possui duas raízes imaginárias dadas por
\[ \lambda = \pm \sqrt{\dfrac{-k}{M}} = \pm i \sqrt{\dfrac{k}{M}} \]
Seja
\[ \omega = \sqrt{\dfrac{k}{M}} \tag{15.6.6} \]
Então, \(\lambda_1 = i\omega\) e \(\lambda_2 = -i\omega\) são as duas raízes. Tanto \(e^{i\omega t}\) quanto \(e^{-i\omega t}\) são soluções da equação diferencial (15.6.3). A solução geral é
\[ A_1 e^{i\omega t} + A_2 e^{-i\omega t} \tag{15.6.7} \]
com constantes arbitrárias reais ou complexas \(A_1\) e \(A_2\). Em geral, esta é uma função complexa do tempo. No entanto, para nosso problema mecânico, apenas uma função real é útil. Portanto, decompomos (15.6.7) em uma parte real e uma parte imaginária e escolhemos as constantes \(A_1\) e \(A_2\) de tal forma que a parte imaginária desapareça. A parte real será então nossa função \(x = x(t)\).
Com base na fórmula (15.4.11), podemos substituir \(e^{i\omega t}\) e \(e^{-i\omega t}\) por combinações lineares de \(\cos(\omega t)\) e \(\sin(\omega t)\). Assim, a parte real de (15.6.7) assume a forma
\[ x = a \cos(\omega t) + b \sin(\omega t) \tag{15.6.8} \]
sendo que \(a\) e \(b\) são constantes reais arbitrárias. É conveniente reescrever essa função na forma
\[ \boxed{x = c \cos\left(\omega \left(t - t_0\right)\right)} \quad (c > 0) \tag{15.6.9} \]
uma função que introduzimos na Seção 5.9. Lá, chamamos \(\omega\) de frequência angular, \(c\) de amplitude e \(t_0\) de acrofase. As relações entre as constantes \(a, b, c, \omega\) e \(t_0\) são dadas por
\[ a = c \cos(\omega t_0), \quad b = c \sin(\omega t_0) \tag{15.6.10} \]
Um gráfico de \(x(t)\) é mostrado na Fig. 5.21c. A equação diferencial (15.6.3) leva a uma oscilação senoidal. A inércia do corpo mantém o sistema em movimento e, na ausência de atrito, a amplitude permanece constante. Por essa razão, o sistema oscilatório é chamado de oscilador harmônico não amortecido.
Agora tornaremos nosso sistema oscilatório mais realista adicionando atrito, como a resistência do ar. Para velocidades moderadas, é um fato experimental que a força de atrito é proporcional à velocidade \(v = dx/dt\). Como o atrito reduz a velocidade, a força de atrito tem sinal oposto ao da velocidade \(v\). Assim, a força de atrito é
\[ -fv = - f \dfrac{dx}{dt} \]
sendo que \(f\) é uma constante positiva. Para obter a força total \(F(x)\) atuando sobre o corpo, devemos adicionar esse termo ao lado direito de (15.6.1). Assim,
\[ F(x) = -kx - f \dfrac{dx}{dt} \tag{15.6.11} \]
A partir de (15.6.2) e (15.6.11), derivamos a equação diferencial
\[ \boxed{M \dfrac{d^2 x}{dt^2} + f \dfrac{dx}{dt} + kx = 0} \quad (M > 0, f > 0, k > 0) \tag{15.6.12} \]
Essa equação também é de segunda ordem e linear, com coeficientes constantes. Usando a função de teste (15.6.4), obtemos
\[ M \lambda^2 A e^{\lambda t} + f \lambda A e^{\lambda t} + k A e^{\lambda t} = 0. \]
Essa equação é satisfeita se, e somente se, \(A = 0\) (caso trivial) ou
\[ \boxed{M \lambda^2 + f \lambda + k = 0} \tag{15.6.13} \]
A equação quadrática (15.6.13) é chamada de equação característica. As duas raízes \(\lambda_1\) e \(\lambda_2\) são dadas por
\[ \lambda = \dfrac{1}{2M} \left( - f \pm \sqrt{f^2 - 4Mk} \right) \tag{15.6.14} \]
A natureza da solução depende do discriminante
\[ D = f^2 - 4Mk \]
Se \(D > 0\), ou seja, se a constante de atrito \(f\) for relativamente grande, as raízes (15.6.14) são reais. Como
\[ \sqrt{f^2 - 4Mk} < f \]
as raízes \(\lambda_1\) e \(\lambda_2\) são ambas negativas. Portanto, a solução da equação diferencial é da forma
\[ c_1 e^{|\lambda_1| t} + c_2 e^{|\lambda_2| t} \]
Essa função de \(t\) é monotonamente decrescente, o que significa que o corpo se desloca lentamente em direção ao seu equilíbrio.
Como estamos interessados em oscilações, descartamos o caso \(D \geq 0\) e exigimos, portanto, que
\[ D = f^2 - 4Mk < 0 \tag{15.6.15} \]
Agora, obtemos duas raízes complexas conjugadas:
\[ \lambda = \dfrac{1}{2M} \left( - f \pm i \sqrt{|D|} \right) \]
Denotamos as raízes por \(-\alpha + i\omega\) e \(-\alpha - i\omega\), respectivamente. A constante \(\alpha = f / 2M\) é positiva. Além disso, temos
\[ \omega = \dfrac{1}{2M} \sqrt{|D|} = \dfrac{1}{2M} \sqrt{4Mk - f^2} \tag{15.6.16} \]
Com essa nova notação, a solução geral da equação diferencial
(15.6.12) é
\[
A_1 e^{(-\alpha + i\omega)t} + A_2 e^{(-\alpha - i\omega)t} = e^{-\alpha
t} \left( A_1 e^{i\omega t} + A_2 e^{-i\omega t} \right)
\]
Como precisamos restringir-nos a uma solução real, podemos tratar a expressão entre parênteses da mesma forma que na equação (15.6.7). Assim, o deslocamento \(x(t)\) assume a forma
\[ x(t) = c e^{-\alpha t} \cos(\omega (t - t_0)) \quad (\alpha > 0) \tag{15.6.17} \]
Comparando essa solução com (15.6.9), observamos que \(x(t)\) ainda oscila com frequência angular \(\omega\), mas a amplitude \(c e^{-\alpha t}\) é monotonamente decrescente. Portanto, o movimento se atenua à medida que \(t\) tende ao infinito. Esse sistema mecânico é chamado de oscilador harmônico amortecido.
Aplicando a fórmula (5.9.2), obtemos para o período \(l\) o valor constante
\[
l = \dfrac{2\pi}{\omega}
\]
Um gráfico de \(x(t)\) é mostrado na Fig. 15.10.
Fig. 15.10. Graph of the displacement x(t) in the case of a damped harmonic oscillator. The amplitude is decreasing. The acrophase is denoted by t0 and the period by l.
Existe uma conexão importante entre a equação diferencial (15.6.12) e o sistema de equações diferenciais de primeira ordem que analisamos na Seção 11.7. Para encontrar essa conexão, substituímos \(dx/dt\) por \(v = v(t)\) e \(d^2x/dt^2\) por \(dv/dt\) na equação (15.6.12). Assim, obtemos
\[ \begin{align} \dfrac{dx}{dt} &= v \\ \dfrac{dv}{dt} &= -\dfrac{k}{M} x - \dfrac{f}{M} v. \end{align} \tag{15.6.18} \]
Essas duas equações têm a forma da equação (11.7.1). Basta identificarmos \(y\) com \(v\) e definirmos
\[ a = 0, \quad b = 1, \quad c = -\dfrac{k}{M}, \quad d = -\dfrac{f}{M} \]
Por outro lado, ao eliminarmos \(v\) de (15.6.18), obtemos a equação de segunda ordem (15.6.12). Assim, (15.6.12) e (15.6.18) são declarações equivalentes.
Deixamos para o leitor a demonstração de que a equação característica (11.7.5) se reduz a (15.6.13).
Os osciladores harmônicos amortecidos não ocorrem apenas na mecânica. Para dar ao leitor uma ideia da ampla aplicabilidade desse conceito, mencionamos um oscilador elétrico típico (Fig. 15.11). Uma corrente elétrica \(I = I(t)\) oscila em um circuito aberto que contém um capacitor de capacitância \(C\), um resistor de resistência \(R\) e uma bobina de autoindutância \(L\). Usando as leis da eletricidade, pode-se derivar a seguinte equação diferencial:
\[ L \dfrac{d^2 I}{dt^2} + R \dfrac{dI}{dt} + \dfrac{1}{C} I = 0 \tag{15.6.19} \]
Fig. 15.11. Um oscilador de circuito elétrico. Ele consiste em um capacitor (capacitância C em farads), uma bobina (indutância L em henries) e um resistor (resistência R em ohms) em série. A corrente elétrica I = I(t) é alternada.
Esta equação é formalmente a mesma que a Eq. (15.6.12). Portanto, se uma condição equivalente à (15.6.15) for satisfeita, a solução será da forma (15.6.17). Comparando o oscilador elétrico com o análogo mecânico, vemos que o resistor desempenha o papel de uma força de atrito, a bobina desempenha o papel da massa inercial e o capacitor desempenha o papel da mola elástica.
Uma generalização adicional do modelo do oscilador é possível permitindo a interferência de forças externas. Por forças externas, entendemos forças que não dependem do estado atual do oscilador. Do ponto de vista biológico, uma força externa pode se originar fora de um organismo, mas também pode ser causada pelo próprio organismo e imposta ao oscilador.
Seja \(E(t)\) a força externa. Então, a equação diferencial (15.6.12) se transforma em:
\[ M \dfrac{d^2 x}{dt^2} + f \dfrac{dx}{dt} + kx = E(t) \tag{15.6.20} \]
A nova equação diferencial ainda é linear em \(x\), mas não é mais homogênea. Chamamos a Eq. (15.6.12) de homogênea, pois, para cada solução \(x(t)\), um múltiplo constante \(A x(t)\) também é uma solução. A nova Eq. (15.6.20) perde essa propriedade se \(E(t) \neq 0\) para alguns valores de \(t\) e, portanto, é chamada de não homogênea. Em biologia, chamamos o sistema de oscilação livre quando \(E(t) = 0\) para todos os valores de \(t\) e forçado no caso oposto.
Equações diferenciais não homogêneas são muito mais difíceis de resolver do que suas contrapartes homogêneas. Por essa razão, consideramos apenas um caso especial. Retornamos às oscilações sem amortecimento, omitindo o atrito. Utilizando a fórmula (15.6.6) e substituindo \(k/M\) por \(\omega^2\) na Eq. (15.6.3), obtemos a equação diferencial homogênea:
\[ \dfrac{d^2 x}{dt^2} + \omega^2 x = 0 \]
Além disso, exigimos que a força externa seja harmônica com amplitude \(A\) e frequência angular \(\omega_1\), ou seja,
\[ E(t) = A \cos(\omega_1 t) \tag{15.6.21} \]
Então, obtemos a equação diferencial não homogênea:
\[ \dfrac{d^2 x}{dt^2} + \omega^2 x = r \cos(\omega_1 t) \tag{15.6.22} \]
onde \(r\) representa \(A/M\). Limitamo-nos a apresentar a solução:
\[ x(t) = c \cos(\omega (t - t_0)) + \dfrac{r}{\omega^2 - \omega_1^2} \cos(\omega_1 t) \tag{15.6.23a} \]
se \(\omega_1 \neq \omega\), e
\[ x(t) = c \cos(\omega (t - t_0)) + \dfrac{rt}{2\omega} \sin(\omega t) \tag{15.6.23b} \]
se \(\omega_1 = \omega\). As constantes \(c\) e \(t_0\) são arbitrárias.
O leitor é convidado a verificar a solução substituindo \(x(t)\) na equação (15.6.22).
Os dois casos \(\omega_1 \neq \omega\) e \(\omega_1 = \omega\) devem ser cuidadosamente distinguidos:
A fórmula (15.6.23a) consiste em dois termos cossenoidais com amplitudes e frequências angulares diferentes. Se \(\omega_1\) difere apenas ligeiramente de \(\omega\), então \(\omega^2 - \omega_1^2\) é um valor pequeno, e a amplitude \(r/(\omega^2 - \omega_1^2)\) é grande. Essa amplitude pode se tornar tão grande que o primeiro termo, com frequência \(\omega\), se torna desprezível. Nesse caso, \(x(t)\) representa aproximadamente uma oscilação harmônica não amortecida. Sua frequência angular é \(\omega_1\), ou seja, a frequência da força externa, e não a frequência \(\omega\) do oscilador livre. Dizemos, portanto, que o oscilador está sincronizado ou treinado pela força externa. O processo de sincronização ou treinamento parece ser importante em sistemas biológicos, onde os períodos em sistemas de oscilação livre podem ser ligeiramente alterados por forças harmônicas externas.
Aqui, a frequência da força externa é exatamente igual à frequência
do oscilador livre. A fórmula (15.6.23b) contém dois termos harmônicos
com a mesma frequência angular \(\omega\). No entanto, o segundo termo
\(\dfrac{rt}{2\omega}\) tende ao
infinito à medida que \(t\)
aumenta. Assim, o sistema é um oscilador harmônico com amplitude
crescente. Portanto, o sistema é instável. Esse fenômeno particular é
conhecido como ressonância.
Para mais detalhes, consulte Defares et al. (1973), Klotter (1960),
Milhorn (1966), Pavlidis (1973), Rosen (1967). Aplicações biológicas são
tratadas por Glass e L. Mackey (1975, p. 136ff.), K. S. Cole (1968,
p. 9), J. M. Smith (1968) e Wiener (1965).
Na Seção 11.8, tratamos de um modelo não linear, que também resultou em oscilações. Esse modelo, conhecido como modelo de Lotka-Volterra, teve sua origem na ecologia. É provável que pesquisas futuras sobre ritmos biológicos envolvam modelos não lineares (cf. Goodwin, 1963, p. 34; Pavlidis, 1973; Wever, 1972).