Aula 8: Função de duas ou mais variáveis independentes
MSP4092: Biomatemática: aspectos quantitativos da vida
17 agosto 2025 17:04h
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(animation, warn.conflicts=FALSE))
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(Deriv, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(nlsr, warn.conflicts=FALSE))
suppressMessages(library(plotly, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(rayshader, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))Adicionar
- Modeling Life - The Mathematics of Biological Systems - Garfinkel et
al - 2017.pdf
- linear approximation to the curve f at the point (X0, Y0) (Figure 7.15), p. 379
- Exercicios 7 e 8
- Explora \(f(X,Y)=e^{−X^2−Y^2}\),
p. 423: use linear approximation relacionado com gradiente
- Use this method to classify the critical points of the functions in Exercise 7.7.10 as local maxima, local minima, or saddle points., p. 431
- Therefore, the Hessian matrix is always symmetric. We can therefore apply a theorem from linear algebra that says that a symmetric matrix can have only real eigenvalues. This has the consequence that there can be no spiraling in a gradient vector field: the state point must head straight upward by the steepest path., p. 431
- 5.2. THE SIR MODEL OF INFECTIOUS DISEASE, p. 360:
- Explorations of mathematical models in biology with MATLAB - Shahin - 2014.pdf
- Fundamentals of Applied Statistics - 4e - Gupta & Kapoor.pdf
- Cap. 4: Demand analysis: elastícidade simple e cruzada
- Cap. 9: Vital Statistics
- Differential Equations and Mathematical Biology - Jones & Sleeman - 2003.pdf
Material
- HTML de R Markdown em
RPubs - Prof. José Siqueira: ResearchGate
Conteúdo
- Número real (Capítulo 1)
- Função e relação (Capítulo 3)
- Funções potência, periódica, exponencial e logarítmica (Capítulos 4, 5 e 6)
- Método gráfico (Capítulo 7)
- Série e limite (Capítulo 8)
- Derivada e integral (Capítulos 9 e 10)
- Equação diferencial ordinária (ODE) (Capítulo 11)
8. Função de duas ou mais variáveis independentes (Capítulo 12)
- Probabilidade (Capítulo 13)
- Matriz, vetor e número complexo (Capítulos 14 e 15)
Introdução
Ceteris paribus é o método de analisar o efeito de uma variável independente sobre a variável dependente mantendo todas as outras variáveis independentes constantes.
A derivada parcial é uma maneira de implementar ceteris paribus na análise de modelo matemático.
Gráfico 3D
Fig. Gráfico 3D: 4 pontos.
Fig. Gráfico de z = f(x, y).
Fig. Algumas superfícies em 3D
Curva de nível
As curvas de nível aparecem em diversas aplicações. Por exemplo, na medicina, se a produção de uma substância \(Q(x, y)\) em um processo biológico for determinada por duas variáveis de entrada \(x\) e \(y\) (como a concentração de um reagente e a atividade enzimática), então a curva de nível \(Q(x, y) = C\) é chamada de curva de produto constante ou, mais brevemente, uma isoquanta (do grego iso, que significa “igual”).
Outra aplicação das curvas de nível na biologia envolve o conceito de curvas de indiferenciação. Por exemplo, ao considerar o comportamento de um organismo que está respondendo a dois estímulos ambientais (como temperatura e umidade), podemos associar uma função de resposta \(U(x, y)\) que mede a satisfação ou adaptação do organismo diante das variáveis. A curva de nível \(U(x, y) = C\) dessa função é chamada de curva de indiferença e fornece todas as combinações de \(x\) e \(y\) que resultam no mesmo nível de resposta ou adaptação do organismo.
Fig. Uma curva de nível da superfície z = f(x, y)
Fig. (a) Superfície z = f(x, y) de uma montanha e (b) curvas de nível produzem um mapa topográfico de z = f(x, y).
Fig. Curvas de nível ajudam a visualizar o formato da superfície.
Geometria da derivada parcial
As derivadas parciais de uma função de duas variáveis podem ser interpretadas geometricamente da seguinte forma. Para cada valor fixo \(y_0\), os pontos \((x, y_0, z)\) formam um plano vertical cuja equação é \(y = y_0\). Se \(z = f(x, y)\) e \(y\) é mantido fixo em \(y = y_0\), então os pontos correspondentes \((x, y_0, f(x, y_0))\) formam uma curva no espaço tridimensional, que é a interseção da superfície \(z = f(x, y)\) com o plano \(y = y_0\).
Em cada ponto dessa curva, a derivada parcial é simplesmente a inclinação da linha no plano \(y = y_0\) que é tangente à curva no ponto em questão. Ou seja, é a inclinação da reta tangente na direção \(x\).
De maneira semelhante, se \(x\) for mantido fixo em \(x = x_0\), os pontos correspondentes \((x_0, y, f(x_0, y))\) formam uma curva que é a interseção da superfície \(z = f(x, y)\) com o plano vertical \(x = x_0\).
Em cada ponto dessa curva, a derivada parcial é a inclinação da linha tangente no plano \(x = x_0\). Ou seja, é a inclinação da reta tangente na direção \(y\).
Fig. Interpretação de derivada parcial.
Exemplo: Média geométrica de duas variáveis
Recordamos a fórmula para da média geométrica de dois números \(x\) e \(y\):
\[ z(x,y) = \sqrt{xy} \\ x \ge 0 \\ y \ge 0 \]
Considere \(x\) e \(y\) como variáveis cujos valores podem ser escolhidos independentemente um do outro. Então, a cada par \((x, y)\) está associado exclusivamente um número \(z\), a média geométrica. No Capítulo 3, chamamos essa associação de função.
Dizemos que \(z\) é uma função do par \((x, y)\), ou o par \((x, y)\) é mapeado em \(z\). Também é comum chamar \(z\) de função de duas variáveis \(x\) e \(y\).
Os números \(x\) e \(y\) são conhecidos como variáveis independentes, enquanto \(z\) é chamado de variável dependente. O domínio \(D\) da função é o conjunto de todos os pares \((x, y)\) com \(x \ge 0\) e \(y \ge 0\).
O domínio é:
\[ D= \{(x, y)|x\ge0, y\ge0\} \]
A imagem \(R\) da função é o conjunto de todos os números \(z \ge0\) ou
\[ R = \{z|\,z\ge 0\} \]
Para um gráfico de uma função de duas variáveis independentes, podemos usar um sistema de coordenadas tridimensionais com eixos \(x\), \(y\), \(z\) perpendiculares aos pares. Um par \((x, y)\) é representado por um ponto no plano \(xy\). A cada um desses pontos está associada uma coordenada \(z\) plotada em uma linha que passa pelo ponto e plotada perpendicularmente ao plano \(xy\). Esta coordenada determina um ponto. Denotamos esse ponto simplesmente por \((x, y, z)\).
Se \(z\) é uma função de \((x, y)\), os pontos \((x, y, z)\) formam uma superfície que pode ser contínua ou não. Normalmente, é desenhada uma vista em perspectiva do sistema de coordenadas e da superfície. A Fig. 12.1 representa a superfície definida pela função \(z\).
Fig. 12.1. Gráfico da média geométrica z = sqrt(x y) em um sistema de coordenadas retangulares.
Nem sempre é conveniente plotar ou usar uma vista em perspectiva. Para funções contínuas, um gráfico de isolinhas pode ser mais fácil de preparar e interpretar. Isolinhas são linhas ao longo das quais o valor funcional permanece constante, da mesma forma que linhas de contorno em um mapa geográfico conectam pontos de elevação constante. Um gráfico de isolinhas da função \(z\) é apresentado na Fig. 12.2.
Fig. 12.2. Isolinhas representando a função z = sqrt(x y).
z <- function(x, y) {
sqrt(x*y)
}
x <- seq(0, 5, length.out = 400)
y <- seq(0, 5, length.out = 400)
grid <- expand.grid(x = x, y = y)
grid$z <- z(grid$x, grid$y)
graf <- ggplot2::ggplot(grid, ggplot2::aes(x, y, z = z)) +
ggplot2::geom_contour(ggplot2::aes(colour = ggplot2::after_stat(level)),
bins = 10) +
ggplot2::labs(title = "Isolinhas de sqrt(x y)",
x = "x",
y = "y",
color = "z") +
ggplot2::theme_minimal()
print(graf)
# Definindo a nova função com sqrt(x * y)
sqrt_xy <- function(a, b) {
sqrt(a * b)
}
# Gerando os valores de a e b
a_values <- seq(0, 5, length.out = 100)
b_values <- seq(0, 5, length.out = 100)
# Criando a grade de combinações de a e b
grid <- base::expand.grid(a = a_values, b = b_values)
# Calculando os valores da função sqrt(a * b)
grid$z <- base::with(grid, sqrt_xy(a, b))
# Transformando os dados para o formato longo
grid_melt <- reshape2::melt(grid, id.vars = c("a", "b"), measure.vars = "z")
# Criando o gráfico de contorno
ggplot_contour <- ggplot2::ggplot(grid_melt, ggplot2::aes(x = a, y = b, z = value)) +
ggplot2::geom_tile(ggplot2::aes(fill = value)) +
ggplot2::geom_contour(color = "black") +
ggplot2::scale_fill_gradientn(colours = grDevices::terrain.colors(10)) +
ggplot2::labs(title = "Gráfico de Contorno da Função sqrt(a * b)", x = "a", y = "b", fill = "C") +
ggplot2::theme_minimal()
# Exibindo o gráfico
print(ggplot_contour)
# Renderizando o gráfico em 3D usando rayshader
rayshader::plot_gg(ggplot_contour, multicore = TRUE, zoom = 0.8, phi = 40, theta = 30)
rayshader::render_snapshot(clear = FALSE)
# Criando o gráfico 3D interativo com plotly
fig <- plotly::plot_ly(
x = ~a_values,
y = ~b_values,
z = ~matrix(grid$z, nrow = 100)
)
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 = "a"),
yaxis = list(title = "b"),
zaxis = list(title = "C"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = -0.64)
)
),
title = "Superfície 3D da Função sqrt(a * b)"
)
fig
x <- seq(0, 5, length.out = 30)
y <- seq(0, 5, length.out = 30)
z <- outer(x, y, function(x, y) sqrt(x * y))
persp(x, y, z, main = "Gráfico de sqrt(x*y)",
xlab = "x", ylab = "y", zlab = "sqrt(x*y)",
theta = 30, phi = 20, expand = 1, col = "lightblue")
persp(x, y, z, main = "Gráfico de sqrt(x*y)",
xlab = "x", ylab = "y", zlab = "sqrt(x*y)",
theta = -30, phi = 20, expand = 1, col = "lightblue")
persp(x, y, z, main = "Gráfico de sqrt(x*y)",
xlab = "x", ylab = "y", zlab = "sqrt(x*y)",
theta = 120, phi = 0, expand = 1, col = "lightblue")
f <- function(ab) {
a <- ab[1]; b <- ab[2]
sqrt(a*b)
}
pracma::ezcontour(f,
xlim = c(0, 5),
ylim = c(0, 5))
Derivada parcial
O símbolo \(\partial\) (conhecido como \(d\) torto ou del) é utilizado para representar derivadas parciais.
Ele é usado quando a função tem mais de uma variável independente, indicando que a derivada é tomada em relação a uma das variáveis, enquanto as outras são mantidas constantes.
Por que usar o \(\partial\) em vez do \(d\)?
O \(d\) comum é usado para derivadas totais, ou seja, quando a função depende de uma única variável independente:
\[ \dfrac{d f}{d x} \]
O \(\partial\) é usado para derivadas parciais, quando a função depende de várias variáveis:
\[ \dfrac{\partial f}{\partial x} \]
Isso indica que estamos derivando \(f\) em relação a \(x\), mantendo as outras variáveis constantes (ceteris paribus).
Exemplo
Seja uma função de duas variáveis:
\[ f(x, y) = x^2 + xy + y^2 \]
As derivadas parciais são:
- Em relação a \(x\):
\[ \dfrac{\partial f}{\partial x} = 2x + y \]
- Em relação a \(y\):
\[ \dfrac{\partial f}{\partial y} = x + 2y \]
Notação de derivada parcial
\[ \begin{align} \dfrac{\partial f(x,y)}{\partial x}&=f_{x}(x,y)=f_{x}\\ \dfrac{\partial f(x,y)}{\partial y}&=f_{y}(x,y)=f_{y}\\ \dfrac{\partial^2 f(x,y)}{\partial x^2}&=f_{xx}(x,y)=f_{xx}\\ \dfrac{\partial^2 f(x,y)}{\partial y^2}&=f_{yy}(x,y)=f_{yy}\\ \dfrac{\partial^2 f(x,y)}{\partial x \partial y}&=\dfrac{\partial^2 f(x,y)}{\partial y \partial x}=f_{xy}(x,y)=f_{xy} \end{align} \]
Exemplo 12.2.1: Média geométrica
Consideramos novamente a média geométrica \(z(x,y)=\sqrt{xy}\).
Mantendo \(y\) constante, \(z\) é proporcional a \(\sqrt{x}\).
Portanto, a derivada parcial resulta em:
\[ \begin{align} \dfrac{\partial z}{\partial x}&=\left(\sqrt{xy}\right)^{\prime} \\ &=\left(x^{1/2}y^{1/2}\right)^{\prime} \\ &=y^{1/2}\left(x^{1/2}\right)^{\prime} \\ &=y^{1/2}\dfrac{1}{2}x^{-1/2} \\ &=\dfrac{1}{2}x^{-1/2}y^{1/2} \\ \dfrac{\partial z}{\partial x} &= \dfrac{1}{2}\sqrt{\dfrac{y}{x}} \end{align} \]
Da mesma forma, quando \(x\) é mantida constante,
\[ \dfrac{\partial z}{\partial y} = \dfrac{1}{2}\sqrt{\dfrac{x}{y}} \]
Exemplo 12.2.2: Corrente elétrica
Consideramos a corrente elétrica \(I\) que flui através de um resistor que é mantido sob temperatura constante.
De acordo com a lei de Ohm,
\[ I=\dfrac{V}{R} \]
sendo que \(V\) denota a tensão e \(R\) a resistência. A corrente \(I\) é uma função de duas variáveis independentes \(V\) e \(R\). Se a tensão \(V\) mudar enquanto \(R\) permanecer constante, obtemos a taxa de variação de:
\[ I_V=\dfrac{\partial I}{\partial V}=\dfrac{1}{R}=R^{-1} \]
Se, no entanto, a resistência \(R\) mudar enquanto \(V\) permanecer constante, a taxa de mudança torna-se
\[ I_R=\dfrac{\partial I}{\partial R}=-\dfrac{V}{R^2}=-VR^{-2} \]
Uma maior diferenciação leva a:
\[ \begin{align} I_{VV}&=\dfrac{\partial^2 I}{\partial^2 V}=0 \\ I_{RR}&=\dfrac{\partial^2 I}{\partial^2 R}=2VR^{-3}=2\dfrac{V}{R^3} \\ I_{VR}&=\dfrac{\partial}{\partial R}\left(\dfrac{\partial I}{\partial V}\right)=\dfrac{\partial}{\partial V}\left(\dfrac{\partial I}{\partial R}\right)=-R^{-2}=-\dfrac{1}{R^2} \end{align} \]
Para obter a derivada mista primeiro diferenciamos em relação a \(V\) e depois em relação a \(R\).
Exemplo: Bens substitutos e complementares
Os bens podem ser classificadas como substitutos ou complementares no mercado.
Os bens substitutos são, por exemplo, manteiga e margarina. Já os bens complementares são, por exemplo, café e açúcar.
Dois bens são substitutos se, ao aumentar o preço de um bem, a demanda pela outra aumenta.
Dois bens são complementares se, ao aumentar o preço de um bem, a demanda pela outra diminui.
Para formalizar e distinguir matematicamente bens complementares e substitutos, é necessário usar a elasticidade cruzada.
A elasticidade instantânea cruzada utiliza a derivada parcial cruzada. Um exemplo é necessário para ilustrar o problema.
As duas funções de demanda, uma para cada bem, são funções dos preços dos dois bens:
\[ \begin{align} d^A&=1000-2p^A+10p^B\\ d^B&=2000+5p^A-20p^B \end{align} \]
sendo \(d^A>0\) e \(d^B>0\) as quantidades semanais em quilograma de manteiga (A) e margarina (B) em um hipermercado, respectivamente, e \(p^A>0\) e \(p^B>0\) são preços unitários da manteiga e da margarina, respectivamente.
Verifique se as funções de demanda são economicamente válidas e se os bens são substitutos.
Solução:
As elasticidades instantâneas parciais são dadas por:
\[ \begin{align} E^A_A&=p^A\dfrac{d^A_A}{d^A}\\ E^B_B&=p^B\dfrac{d^B_B}{d^B} \end{align} \]
O sinal das elasticidades instantâneas parciais dependem diretamente das derivadas parciais das funções de demanda.
As derivadas parciais são dadas por:
\[ \begin{align} d^A_A&=-2\\ d^B_B&=-20 \end{align} \]
Como essas derivadas parciais são estritamente negativas, as funções de demanda são economicamente plausíveis (aumento de preço do bem causa diminuição de demanda deste bem) e, portanto, as elasticidades instantâneas parciais são estritamente negativas:
\[ \begin{align} E^A_A&=-\dfrac{2p^A}{d^A}<0\\ E^B_B&=-\dfrac{20p^B}{d^B}<0 \end{align} \]
As derivadas parciais cruzadas são:
\[ \begin{align} d^A_B&=10\\ d^B_A&=5 \end{align} \]
Observe que as duas derivadas parciais cruzadas são estritamente positivas.
As elasticidades instantâneas cruzadas da demanda pelos bens A e B são expressas, respectivamente, da seguinte maneira:
\[ \begin{align} E^A_B&=p^B\dfrac{d^A_B}{d^A}\\ E^B_A&=p^A\dfrac{d^B_A}{d^B} \end{align} \]
O sinal da elasticidade instantânea cruzada depende diretamente do sinal da derivada parcial cruzada.
Como essas derivadas parciais cruzadas são estritamente positivas, as funções de demanda são economicamente plausíveis (aumento de preço do bem causa aumento de demanda pelo outro bem) e, portanto, as elasticidades instantâneas parciais cruzadas são estritamente positivas:
\[ \begin{align} E^A_B&=\dfrac{10p^B}{d^A}>0\\ E^B_A&=\dfrac{5p^A}{d^B}>0 \end{align} \]
Dois bens são substitutos se as elasticidades instantâneas cruzadas são estritamente positivas, ou seja:
\[ \begin{align} E^A_B&>0\\ E^B_A&>0 \end{align} \]
Por outro lado, dois bens são complementares se as elasticidades instantâneas cruzadas são estritamente negativas.
\[ \begin{align} E^A_B&<0\\ E^B_A&<0 \end{align} \]
Conclui-se que manteiga e margarina são bens substitutos.
Aproximação de Taylor de primeira ordem
A aproximação linear de Taylor para uma função de duas variáveis \(f(x, y)\) é uma extensão da aproximação de Taylor de uma variável. Ela fornece uma maneira de aproximar o valor da função perto de um ponto conhecido \((a, b)\).
A aproximação linear é dada por:
\[ f(x, y) \approx f(a, b) + \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) \]
Aqui:
- \(f(a, b)\) é o valor da função no ponto de expansão \((a, b)\).
- \(\dfrac{\partial f}{\partial x}(a, b)\) e \(\dfrac{\partial f}{\partial y}(a, b)\) são as derivadas parciais da função em relação a \(x\) e \(y\) avaliadas no ponto \((a, b)\).
- \((x - a)\) e \((y - b)\) representam os deslocamentos a partir do ponto \((a, b)\).
A aproximação linear corresponde à equação do plano tangente à superfície dada por \(f(x, y)\) no ponto \((a, b)\). Isso significa que, para valores de \(x\) e \(y\) próximos a \((a, b)\), o valor de \(f(x, y)\) pode ser aproximado pelo plano tangente.
Esta aproximação é especialmente útil quando a função é complexa e precisamos de um valor aproximado rapidamente para cálculos próximos a um ponto específico.
Exemplo: Plano Tangente à Superfície
Para visualizar o plano tangente junto à superfície, siga estas etapas:
Definir a Função: Escolha uma função de duas variáveis, por exemplo, \(f(x, y) = x^2 + y^2\).
Calcular as Derivadas Parciais:
- \(\dfrac{\partial f}{\partial x} = 2x\)
- \(\dfrac{\partial f}{\partial y} = 2y\)
- Plano Tangente no Ponto (a, b): Substitua as derivadas no ponto escolhido, por exemplo, \((1, 1)\):
\[ z = 1^2 + 1^2 + 2 \cdot 1 \cdot (x - 1) + 2 \cdot 1 \cdot (y - 1) \]
Simplificando:
\[ \begin{align} z &= 2 + 2(x - 1) + 2(y - 1) z &= 2x + 2y - 2 \end{align} \]
# Definir a função original
f <- function(x, y) {
x^2 + y^2
}
# Parâmetros do ponto tangente
a <- 1
b <- 1
# Derivadas parciais
df_dx <- function(x, y) { 2 * x }
df_dy <- function(x, y) { 2 * y }
# Plano tangente
tangent_plane <- function(x, y) {
f(a, b) + df_dx(a, b) * (x - a) + df_dy(a, b) * (y - b)
}
# Criar uma grade de pontos
x_vals <- seq(0, 2, length.out = 30)
y_vals <- seq(0, 2, length.out = 30)
# Criar uma malha
grid <- expand.grid(x = x_vals, y = y_vals)
# Calcular os valores da função e do plano tangente
grid$z <- mapply(f, grid$x, grid$y)
grid$tangent <- mapply(tangent_plane, grid$x, grid$y)
# Gráfico 3D com plotly (sem piping)
fig <- plotly::plot_ly()
# Superfície da função - cor azul
fig <- plotly::add_surface(fig, x = ~x_vals, y = ~y_vals, z = matrix(grid$z, nrow = 30),
showscale = FALSE, name = "Função f(x, y)",
surfacecolor = matrix(grid$z, nrow = 30), colorscale = "Blues")
# Plano tangente - cor laranja
fig <- plotly::add_surface(fig, x = ~x_vals, y = ~y_vals, z = matrix(grid$tangent, nrow = 30),
showscale = FALSE, name = "Plano Tangente",
surfacecolor = matrix(grid$tangent, nrow = 30), colorscale = "Oranges")
# Layout do gráfico
fig <- plotly::layout(fig,
title = "Superfície da Função e Plano Tangente",
scene = list(xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "f(x, y) = x² + y²")))
figExemplo:
Para obter a aproximação de Taylor de primeira ordem da função \(\exp(x^2 + y^2)\) em torno do ponto \((a, b)\), utilizamos a seguinte fórmula da série de Taylor de primeira ordem:
\[ f(x, y) \approx f(a, b) + \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) \]
Vamos calcular manualmente a aproximação de Taylor de primeira ordem da função \(\exp(x^2 + y^2)\) em torno do ponto \((a, b)\).
Primeiro, calculamos as derivadas parciais:
\[ \dfrac{\partial}{\partial x} [\exp(x^2 + y^2)] = 2x \cdot \exp(x^2 + y^2) \]
\[ \dfrac{\partial}{\partial y} [\exp(x^2 + y^2)] = 2y \cdot \exp(x^2 + y^2) \]
Substituindo \(x = a\) e \(y = b\), temos:
\[ \dfrac{\partial}{\partial x} [\exp(a^2 + b^2)] = 2a \cdot \exp(a^2 + b^2) \]
\[ \dfrac{\partial}{\partial y} [\exp(a^2 + b^2)] = 2b \cdot \exp(a^2 + b^2) \]
A aproximação de Taylor de primeira ordem é dada por:
\[ \exp(x^2 + y^2) \approx \exp(a^2 + b^2) + 2a \cdot \exp(a^2 + b^2) \cdot (x - a) + 2b \cdot \exp(a^2 + b^2) \cdot (y - b) \]
Derivada total: Regra da cadeia generalizada
A derivada total é uma ferramenta matemática poderosa utilizada para analisar a variação de uma função que depende de múltiplas variáveis em relação a uma única variável. Enquanto a derivada parcial mede a taxa de variação de uma função em relação a uma variável específica, mantendo as outras constantes, a derivada total leva em consideração a variação de todas as variáveis independentes simultaneamente.
Quando uma função \(f\) depende indiretamente de uma variável \(t\) através de outras variáveis (\(x\) e \(y\)), temos:
\[ x = x(t) \quad \text{e} \quad y = y(t) \]
A regra da cadeia nos diz que a derivada total de \(f\) com relação a \(t\) é:
\[ \dfrac{df}{dt} = \dfrac{\partial f}{\partial x} \cdot \dfrac{dx}{dt} + \dfrac{\partial f}{\partial y} \cdot \dfrac{dy}{dt} \]
Essa expressão considera que tanto \(x\) quanto \(y\) variam em função de \(t\).
A derivada total é uma aplicação direta da regra da cadeia quando as variáveis independentes dependem de uma outra variável comum. A expressão da derivada total é:
\[ df = \dfrac{\partial f}{\partial x} \, dx + \dfrac{\partial f}{\partial y} \, dy \]
Quando associamos \(dx\) e \(dy\) às variações induzidas por \(t\), temos:
\[ dx = \dfrac{dx}{dt} \, dt \quad \text{e} \quad dy = \dfrac{dy}{dt} \, dt \]
Substituindo na expressão da derivada total:
\[ df = \left[\dfrac{\partial f}{\partial x} \cdot \dfrac{dx}{dt} + \dfrac{\partial f}{\partial y} \cdot \dfrac{dy}{dt}\right] \, dt \]
Isso mostra que a derivada total incorpora a regra da cadeia quando as variáveis dependem de um parâmetro comum.
Exemplo: Regra da cadeia para computar taxa de demanda
Diná gerencia uma loja de produtos naturais que vende dois tipos de água vitaminada: marca A e marca B. Os registros de vendas indicam que, se a marca A for vendida por \(x\) dólares por garrafa e a marca B por \(y\) dólares por garrafa, a demanda pela marca A será dada por:
\[ Q(x, y) = 300 - 20x^2 + 30y \quad \text{garrafas por mês} \]
Ela estima que, daqui a \(t\) meses, o preço da marca A será:
\[ x(t) = 2 + 0.05t \quad \text{dólares por garrafa} \]
E o preço da marca B será:
\[ y(t) = 2 + 0.1\sqrt{t} \quad \text{dólares por garrafa} \]
Pergunta:
A que taxa Diná deve esperar que a demanda pela marca A esteja mudando em relação ao tempo daqui a 4 meses? Ela deve esperar que a demanda esteja aumentando ou diminuindo nesse momento?
Solução:
O objetivo de Diná é encontrar \(\dfrac{dQ}{dt}\) quando \(t = 4\). Usando a regra da cadeia, temos:
\[ \dfrac{dQ}{dt} = \dfrac{\partial Q}{\partial x} \dfrac{dx}{dt} + \dfrac{\partial Q}{\partial y} \dfrac{dy}{dt} \]
Substituindo as derivadas parciais e as taxas de variação:
\[ \dfrac{dQ}{dt} = -40x\times 0.05 + 30\times(0.05t^{-1/2}) \]
Quando \(t = 4\):
\[ x(4) = 2 + 0.05\times 4 = 2.2 \]
Portanto:
\[ \dfrac{dQ}{dt} = -40(2.2)(0.05) + 30(0.05)(0.5) = -3.65 \]
Conclusão:
Daqui a 4 meses, a demanda mensal pela marca A estará diminuindo a uma taxa de 3,65 garrafas por mês.
Exemplo: Crescimento Populacional com Disponibilidade de Alimento e Nível de Predação
Vamos estudar uma população de coelhos \(P\) que depende da quantidade de alimento disponível \(A\) e nível de predação \(R\). A relação pode ser expressa como uma função:
\[ P = f(A, R) \]
Identificação das Variáveis:
- \(P\): População de coelhos (variável dependente).
- \(A\): Quantidade de alimento disponível (variável independente).
- \(R\): Nível de predação (variável independente).
Função da População:
- Suponhamos que a função da população seja dada por:
\[ P(A, R) = \dfrac{A^2}{R} \]
Derivadas Parciais:
- \(\dfrac{\partial P}{\partial A}\): A taxa de variação da população em relação à quantidade de alimento.
\[ \dfrac{\partial P}{\partial A} = \dfrac{2A}{R} \]
- \(\dfrac{\partial P}{\partial R}\): A taxa de variação da população em relação ao nível de predação.
\[ \dfrac{\partial P}{\partial R} = -\dfrac{A^2}{R^2} \]
Taxas de Variação das Variáveis:
- Suponha que a quantidade de alimento disponível varia de forma senoidal com o tempo:
\[ A = \sin(t) \]
Então, a taxa de variação de $A$ em relação ao tempo é:
\[ \dfrac{dA}{dt} = \cos(t) \]
- E o nível de predação varia como a raiz quadrada do tempo:
\[ R = \sqrt{t} \]
Então, a taxa de variação de $R$ em relação ao tempo é:
\[ \dfrac{dR}{dt} = \dfrac{1}{2\sqrt{t}} \]
- Derivada Total:
A derivada total de \(P\) em relação ao tempo \(t\) é dada por:
\[ \dfrac{dP}{dt} = \dfrac{\partial P}{\partial A} \dfrac{dA}{dt} + \dfrac{\partial P}{\partial R} \dfrac{dR}{dt} \]
Substituindo os valores das derivadas parciais e das taxas de variação:
- \(\dfrac{\partial P}{\partial A} = \dfrac{2A}{R} = \dfrac{2 \sin(t)}{\sqrt{t}}\)
- \(\dfrac{\partial P}{\partial R} = -\dfrac{A^2}{R^2} = -\dfrac{\sin^2(t)}{t}\)
Então:
\[ \dfrac{dP}{dt} = \left( \dfrac{2 \sin(t)}{\sqrt{t}} \right) \cos(t) + \left( -\dfrac{\sin^2(t)}{t} \right) \dfrac{1}{2\sqrt{t}} \]
Portanto, a derivada total da população de coelhos em relação ao tempo \(t\) é:
\[ \dfrac{dP}{dt} = \dfrac{\sin(t)}{\sqrt{t}} \left( 2 \cos(t) - \dfrac{\sin(t)}{2t} \right) \]
Aproximação incremental da derivada total
Em particular, se \(f(x)\), então:
\[ \Delta f \approx \dfrac{df}{dx} \Delta x \]
onde \(\Delta x\) é uma pequena mudança na variável \(x\) e \(\Delta f\) é a mudança correspondente em \(f\).
Aqui está a fórmula análoga de aproximação para funções de duas variáveis, baseada na regra da cadeia para derivadas parciais.
Suponha que \(f\) seja uma função de \(x\) e \(y\). Se \(\Delta x\) denota uma pequena mudança em \(x\) e \(\Delta y\) uma pequena mudança em \(y\), a mudança correspondente em \(f(x,y)\) é aproximada por:
\[ \Delta f \approx \dfrac{\partial f}{\partial x} \Delta x + \dfrac{\partial f}{\partial y} \Delta y \]
Relação com a Aproximação de Taylor e Regra da Cadeia
A fórmula da aproximação de Taylor de primeira ordem para funções de duas variáveis é:
\[ \begin{align} f(x, y) &\approx f(a, b) + \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) \\ f(x, y) - f(a, b) &\approx + \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) \\ \end{align} \]
Podemos reorganizar essa expressão para destacar a variação na função:
\[ \Delta f \approx \dfrac{\partial f}{\partial x} \cdot \Delta x + \dfrac{\partial f}{\partial y} \cdot \Delta y \]
Essa expressão é precisamente a fórmula da derivada total (\(df\)), onde:
- \(\Delta f\) é a variação aproximada da função.
- \(\Delta x\) e \(\Delta y\) correspondem aos pequenos incrementos correspondentes a \(dx\) e \(dy\).
A derivada total fornece a taxa de variação instantânea da função considerando todas as variáveis ao mesmo tempo. A aproximação de Taylor, por sua vez, usa essa informação para estimar o valor da função em torno de um ponto conhecido.
Assim, a aproximação de Taylor de primeira ordem pode ser vista como uma aplicação da derivada total, fornecendo uma forma linear de prever mudanças na função devido a pequenas perturbações nas variáveis.
A relação entre a derivada total, a aproximação de Taylor e a regra da cadeia surge naturalmente quando lidamos com funções compostas e funções de várias variáveis. Vamos explorar essa conexão.
Aproximação de Taylor e Regra da Cadeia
Na aproximação de Taylor de primeira ordem, temos:
\[ \Delta f \approx \dfrac{\partial f}{\partial x} \cdot \Delta x + \dfrac{\partial f}{\partial y} \cdot \Delta y \]
Se \(\Delta x\) e \(\Delta y\) resultam de pequenas variações na variável \(t\), a aproximação se torna:
\[ \Delta f \approx \dfrac{\partial f}{\partial x} \cdot \dfrac{dx}{dt} \cdot \Delta t + \dfrac{\partial f}{\partial y} \cdot \dfrac{dy}{dt} \cdot \Delta t \]
Ou seja, a aproximação de Taylor, quando aplicada a variáveis dependentes de outra, também utiliza a regra da cadeia para determinar a variação da função.
A regra da cadeia é essencial para calcular a derivada total quando as variáveis são funções de outra variável. A aproximação de Taylor de primeira ordem aproveita essa relação para estimar a variação da função em torno de um ponto, especialmente quando as variáveis dependem de um parâmetro comum.
Exemplo: Aproximação Incremental da Produção
Em uma determinada fábrica, a produção diária é dada por:
\[ Q = 60K^{1/2}L^{1/3} \]
onde:
- \(K\) representa o investimento de capital (em milhares de dólares).
- \(L\) representa a força de trabalho (em horas de trabalho).
Atualmente:
- O investimento de capital é de $900.000 (900 unidades).
- A força de trabalho é de 1.000 horas por dia.
Problema:
Estime a mudança na produção se o investimento de capital aumentar em $1.000 (1 unidade) e a força de trabalho aumentar em 2 horas.
Solução:
Aplicando a fórmula de aproximação incremental com:
- \(K = 900\), \(L = 1.000\)
- \(\Delta K = 1\) ($1.000)
- \(\Delta L = 2\)
\[ \Delta Q \approx \dfrac{\partial Q}{\partial K} \Delta K + \dfrac{\partial Q}{\partial L} \Delta L \]
Calculando as derivadas parciais:
\[ \dfrac{\partial Q}{\partial K} = 30K^{-1/2}L^{1/3} \quad \text{e} \quad \dfrac{\partial Q}{\partial L} = 20K^{1/2}L^{-2/3} \]
Substituindo os valores:
\[ \begin{align} \Delta Q &\approx 30(1/30)(10)(1) + 20(30)(1/100)(2)\\ \Delta Q &\approx 10 + 12 = 22 \, \text{unidades} \end{align} \]
Conclusão:
A produção aumentará em aproximadamente 22 unidades.
Exemplo: Aproximação Incremental da Área Superficial do Corpo Humano
Sabemos que a área superficial do corpo humano pode ser medida pela fórmula empírica:
\[ S(W, H) = 0.0072\times W^{0.425}H^{0.725} \]
onde:
- \(W\) (kg) é o peso da pessoa.
- \(H\) (cm) é a altura da pessoa.
Atualmente, uma criança tem massa corportal 34 kg e mede 120 cm de estatura.
Questão a)
Calcule as derivadas parciais \(S_W(34, 120)\) e \(S_H(34, 120)\) e interprete cada uma como uma taxa de variação.
Questão b)
Estime a mudança na área superficial se a massa corporal da criança aumentar em 1 kg e a estatura aumentar em 3 cm.
# Definindo a função
f <- function(x, y) { 0.0072*x^(0.425)*y^(0.725) }
# Derivadas parciais usando o pacote Deriv
df_dx <- Deriv::Deriv(f, "x")
df_dy <- Deriv::Deriv(f, "y")
# Gradiente da função
gradiente <- Deriv::Deriv(f, c("x", "y"))
print(gradiente)function (x, y)
c(x = 0.00306 * (y^0.725/x^0.575), y = 0.00522 * (x^0.425/y^0.275))# Valores iniciais
x0 <- 34
y0 <- 120
Dx <- 1
Dy <- 2
D <- c(Dx, Dy)
# Gradiente no ponto (x0, y0)
grad_val <- do.call(gradiente, list(x = x0, y = y0))
incremental <- crossprod(grad_val, D)
cat("Gradiente: ", grad_val, "\n")Gradiente: 0.01295769 0.006262882 f(x0, y0) = 1.036615 Aproximação Incremental = 0.02548345 f(x0, y0) + Aproximação Incremental = 1.062098 f(x0 + Dx, y0 + Dy) = 1.062117 Exemplo: Lata cilíndrica com tampa
Uma lata cilíndrica com tampa tem raio \(r\) e altura \(h\). As fórmulas do volume, superfície e custo total são:
- \(v(r,h) = \pi r^2 h\)
- \(s(h,r) = 2\pi r^2 + 2\pi rh\)
- \(c(h, r) = 0.3 \times s(h,r) + 0.2 \times v(h,r)\)
A lata tem raio de 3 cm e altura de 12 cm.
Questão a)
Calcule as derivadas parciais do volume, superfície e custo total em relação ao raio e altura e interprete cada uma como uma taxa de variação.
Questão b)
Estime a mudança no volume, superfície e custo total se o raio da lata diminuir em 1 cm e a altura aumentar em 2 cm.
# Definindo a função
f <- function(x, y) { pi*y*x^2 }
# Derivadas parciais usando o pacote Deriv
df_dx <- Deriv::Deriv(f, "x")
df_dy <- Deriv::Deriv(f, "y")
# Gradiente da função
gradiente <- Deriv::Deriv(f, c("x", "y"))
print(gradiente)function (x, y)
c(x = 2 * (pi * x * y), y = pi * x^2)# Valores iniciais
x0 <- 3
y0 <- 12
Dx <- -1
Dy <- 2
D <- c(Dx, Dy)
# Gradiente no ponto (x0, y0)
grad_val <- do.call(gradiente, list(x = x0, y = y0))
incremental <- crossprod(grad_val, D)
cat("Gradiente: ", grad_val, "\n")Gradiente: 226.1947 28.27433 f(x0, y0) = 339.292 Aproximação Incremental = -169.646 f(x0, y0) + Aproximação Incremental = 169.646 f(x0 + Dx, y0 + Dy) = 175.9292 # Definindo a função
f <- function(x, y) { 2*pi*x^2 + 2*pi*x*y }
# Derivadas parciais usando o pacote Deriv
df_dx <- Deriv::Deriv(f, "x")
df_dy <- Deriv::Deriv(f, "y")
# Gradiente da função
gradiente <- Deriv::Deriv(f, c("x", "y"))
print(gradiente)function (x, y)
c(x = pi * (2 * y + 4 * x), y = 2 * (pi * x))# Gradiente no ponto (x0, y0)
grad_val <- do.call(gradiente, list(x = x0, y = y0))
incremental <- crossprod(grad_val, D)
cat("Gradiente: ", grad_val, "\n")Gradiente: 113.0973 18.84956 f(x0, y0) = 282.7433 Aproximação Incremental = -75.39822 f(x0, y0) + Aproximação Incremental = 207.3451 f(x0 + Dx, y0 + Dy) = 201.0619 # Definindo a função
f <- function(x, y) { 0.3*(pi*y*x^2) + 0.2*(2*pi*x^2 + 2*pi*x*y) }
# Derivadas parciais usando o pacote Deriv
df_dx <- Deriv::Deriv(f, "x")
df_dy <- Deriv::Deriv(f, "y")
# Gradiente da função
gradiente <- Deriv::Deriv(f, c("x", "y"))
print(gradiente)function (x, y)
c(x = pi * (0.2 * (2 * y + 4 * x) + 0.6 * (x * y)), y = pi *
x * (0.3 * x + 0.4))# Gradiente no ponto (x0, y0)
grad_val <- do.call(gradiente, list(x = x0, y = y0))
incremental <- crossprod(grad_val, D)
cat("Gradiente: ", grad_val, "\n")Gradiente: 90.47787 12.25221 f(x0, y0) = 158.3363 Aproximação Incremental = -65.97345 f(x0, y0) + Aproximação Incremental = 92.36282 f(x0 + Dx, y0 + Dy) = 92.99114 Vetor gradiente
O conceito de gradiente é fundamental em cálculo multivariado e é amplamente utilizado em otimização e análise de funções. O gradiente de uma função escalar de várias variáveis é um vetor que aponta na direção de maior taxa de aumento da função. Vamos usar a função \(f(x, y) = \sqrt{x y}\) como exemplo para entender como calcular o gradiente.
A função \(f(x, y) = \sqrt{x y}\) pode ser reescrita como \(f(x, y) = (xy)^{1/2}\).
O gradiente de \(f\) é o vetor das suas derivadas parciais:
\[ \nabla f(x, y) = \left( \dfrac{\partial f}{\partial x}, \dfrac{\partial f}{\partial y} \right) \]
Portanto, o gradiente de \(f(x, y) = \sqrt{xy}\) é:
\[ \nabla f(x, y) = \left( \dfrac{y}{2 \sqrt{xy}}, \dfrac{x}{2 \sqrt{xy}} \right) \]
O gradiente \(\nabla f(x, y)\) aponta na direção de maior aumento de \(f\) e sua magnitude indica a taxa de variação de \(f\) nessa direção.
f <- function(x, y) {
sqrt(x * y)
}
grad_f <- function(x, y) {
c(y / (2 * sqrt(x * y)), x / (2 * sqrt(x * y)))
}
# Gerar uma grade de pontos
x <- seq(0.1, 5, length.out = 10)
y <- seq(0.1, 5, length.out = 10)
grid <- expand.grid(x = x, y = y)
# Calcular os valores da função e do gradiente
grid$z <- mapply(f, grid$x, grid$y)
gradients <- t(mapply(grad_f, grid$x, grid$y))
grid$u <- gradients[, 1]
grid$v <- gradients[, 2]
# Criar o gráfico
g <- ggplot2::ggplot(grid, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z), alpha = 0.3) +
ggplot2::geom_segment(ggplot2::aes(xend = x + u,
yend = y + v),
arrow = ggplot2::arrow(length = ggplot2::unit(0.3, "cm")),
color = "black") +
ggplot2::scale_fill_gradient(low = "white", high = "black") +
ggplot2::labs(title = "Vetor Gradiente de f(x, y) = sqrt(x * y)",
x = "x", y = "y") +
ggplot2::theme_minimal()
print(g)
Na Figura 12.3, o gradiente no ponto de máximo \((6,3)\) é nulo, i.e., \(\nabla f(6, 3)=(0,0)\).
Fig. 12.3. Encontrando o máximo de função.
Matriz Hessiana
A matriz Hessiana foi desenvolvida no século XIX pelo matemático alemão Ludwig Otto Hesse e posteriormente nomeada em sua homenagem. Hesse originalmente usou o termo “determinantes funcionais”.
A matriz Hessiana é uma matriz quadrada e simétrica de segundas derivadas parciais de uma função. Ela é fundamental em cálculo multivariado e otimização para entender a curvatura e a natureza dos pontos críticos da função.
Matriz Hessiana descreve a curvatura local de uma função de várias variáveis.
A Hessiana é às vezes denotada por \(\mathcal{H}\) ou \(\nabla^2\).
\[\mathcal{H}(f(x,y))=\left[\begin{array}{cc} f_{xx} & f_{xy}\\ f_{yx} & f_{yy} \end{array}\right] = \nabla(\nabla f(x,y)) = (J(\nabla f(x,y))^T\]
Exemplo
Para determinar a matriz Hessiana de \(\sqrt{x y}\), calculamos todas as derivadas parciais de segunda ordem.
A função dada é \(f(x, y) = \sqrt{x y}\).
Primeiro, calculamos as derivadas parciais de primeira ordem:
\[ \dfrac{\partial f}{\partial x} = \dfrac{\partial}{\partial x} \sqrt{x y} = \dfrac{y}{2 \sqrt{x y}} \]
\[ \dfrac{\partial f}{\partial y} = \dfrac{\partial}{\partial y} \sqrt{x y} = \dfrac{x}{2 \sqrt{x y}} \]
Agora, calculamos as derivadas parciais de segunda ordem:
- Derivada parcial de \(\dfrac{\partial f}{\partial x}\) em relação a \(x\):
\[ \dfrac{\partial}{\partial x} \left( \dfrac{y}{2 \sqrt{x y}} \right) = -\dfrac{y^2}{4 (x y)^{3/2}} \]
- Derivada parcial de \(\dfrac{\partial f}{\partial x}\) em relação a \(y\):
\[ \dfrac{\partial}{\partial y} \left( \dfrac{y}{2 \sqrt{x y}} \right) = \dfrac{1}{2 \sqrt{x y}} - \frac{1}{4 (x y)^{1/2}} \]
Como \(\dfrac{x y}{(x y)^{3/2}} = \dfrac{1}{(x y)^{1/2}}\), temos:
\[ \dfrac{\partial^2 f}{\partial y \partial x} = \dfrac{1}{4 \sqrt{x y}} \]
- Derivada parcial de \(\dfrac{\partial f}{\partial y}\) em relação a \(x\):
\[ \dfrac{\partial}{\partial x} \left( \dfrac{x}{2 \sqrt{x y}} \right) = \dfrac{1}{2 \sqrt{x y}} - \dfrac{1}{4 (x y)^{1/2}} \]
Como \(\dfrac{x y}{(x y)^{3/2}} = \dfrac{1}{(x y)^{1/2}}\), temos:
\[ \dfrac{\partial^2 f}{\partial x \partial y} = \dfrac{1}{4 \sqrt{x y}} \]
- Derivada parcial de \(\dfrac{\partial f}{\partial y}\) em relação a \(y\):
\[ \dfrac{\partial}{\partial y} \left( \dfrac{x}{2 \sqrt{x y}} \right) = -\dfrac{x^2}{4 (x y)^{3/2}} \]
Portanto, a matriz Hessiana de \(\sqrt{x y}\) é:
\[ \mathcal{H} = \begin{pmatrix} -\dfrac{y^2}{4 (x y)^{3/2}} & \dfrac{xy}{4 (x y)^{3/2}} \\ \dfrac{xy}{4 (x y)^{3/2}} & -\dfrac{x^2}{4 (x y)^{3/2}} \end{pmatrix} \]
Para encontrar a matriz Hessiana de \(\sqrt{x y}\) com \(x = 4\) e \(y = 1\), substituímos esses valores na matriz Hessiana simplificada:
\[ \mathcal{H} = \begin{pmatrix} -\dfrac{1}{32} & \dfrac{1}{8} \\ \dfrac{1}{8} & -\dfrac{1}{2} \end{pmatrix} \]
Portanto, a matriz Hessiana em valores decimais é:
\[ \mathcal{H} = \begin{pmatrix} -0.03125 & 0.125 \\ 0.125 & -0.5 \end{pmatrix} \]
partial derivative sqrt(x y) where x=4, y=1
{D[Sqrt[x y], x], D[Sqrt[x y], y]}
second derivative sqrt(x y) for x and x
second derivative sqrt(x y) for y and y
second derivative sqrt(x y) for x and y
second derivative sqrt(x y) for y and x
{D[Sqrt[x y], {x,2}], D[Sqrt[x y], {y,2}]}
{D[Sqrt[x y], x,y}], D[Sqrt[x y], y,x]}
D[Sqrt[x y], x , y] /. x=4, y=1
Cálculo multivariado em R
# symbolic partial derivatives
f <- expression(sqrt(x*y))
f_x <- calculus::derivative(f=f, var="x", order=1)
f_x[1] "0.5 * (y * (x * y)^-0.5)"[1] "0.5 * (x * (x * y)^-0.5)"# symbolic higher order partial derivatives
f_xx <- calculus::derivative(f=f, var="x", order=2)
f_xx[1] "0.5 * (y * (-0.5 * (y * (x * y)^-1.5)))"[1] "0.5 * (x * (-0.5 * (x * (x * y)^-1.5)))"[1] "0.5 * ((x * y)^-0.5 + y * (-0.5 * (x * (x * y)^-1.5)))"[1] "0.5 * (y * (x * y)^-0.5)" "0.5 * (x * (x * y)^-0.5)" [,1]
[1,] "-(0.5 * (y * ((x * y)^-(0.5 + 1) * (0.5 * y))))"
[2,] "0.5 * ((x * y)^-0.5 - x * ((x * y)^-(0.5 + 1) * (0.5 * y)))"
[,2]
[1,] "0.5 * ((x * y)^-0.5 - y * ((x * y)^-(0.5 + 1) * (0.5 * x)))"
[2,] "-(0.5 * (x * ((x * y)^-(0.5 + 1) * (0.5 * x))))" [1] 0.000000 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751
[9] 2.828427 3.000000 3.162278 [1] 0.000000 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751
[9] 2.828427 3.000000 3.162278
attr(,"gradient")
x y
[1,] Inf NaN
[2,] 1.0000000 0.2500000
[3,] 0.7071068 0.3535534
[4,] 0.5773503 0.4330127
[5,] 0.5000000 0.5000000
[6,] 0.4472136 0.5590170
[7,] 0.4082483 0.6123724
[8,] 0.3779645 0.6614378
[9,] 0.3535534 0.7071068
[10,] 0.3333333 0.7500000
[11,] 0.3162278 0.7905694f_xn <- calculus::derivative(f=f, var=c(x=4, y=1), order=1)
f_xn # calcula primeiro a derivada parcial em x [,1] [,2]
[1,] 0.25 1 [,1] [,2]
[1,] 0.25 1 [,1] [,2]
[1,] -0.03125 -0.5 [,1] [,2]
[1,] -0.03125 -0.5[1] 0.125[1] 0.25 1.00 [,1] [,2]
[1,] -0.03125 0.125
[2,] 0.12500 -0.500Aproximacao de Taylor de segunda ordem
A aproximação de Taylor de segunda ordem para uma função de duas variáveis \(f(x, y)\) é uma expansão mais precisa que a de primeira ordem, pois inclui os termos quadráticos.
A expansão de Taylor de segunda ordem em torno do ponto \((a, b)\) é dada por:
\[ f(x, y) \approx f(a, b) + \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) +\\ +\dfrac{1}{2} \dfrac{\partial^2 f}{\partial x^2}(a, b) \cdot (x - a)^2 + \dfrac{1}{2} \dfrac{\partial^2 f}{\partial y^2}(a, b) \cdot (y - b)^2 +\\ + \frac{\partial^2 f}{\partial x \partial y}(a, b) \cdot (x - a)(y - b) \]
Componentes da Fórmula:
Termo constante: \(f(a, b)\)
Termos lineares:
- \(\dfrac{\partial f}{\partial x}(a, b) \cdot (x - a)\)
- \(\dfrac{\partial f}{\partial y}(a, b) \cdot (y - b)\)
Termos quadráticos:
- \(\dfrac{1}{2} \dfrac{\partial^2 f}{\partial x^2}(a, b) \cdot (x - a)^2\)
- \(\dfrac{1}{2} \dfrac{\partial^2 f}{\partial y^2}(a, b) \cdot (y - b)^2\)
- \(\dfrac{\partial^2 f}{\partial x \partial y}(a, b) \cdot (x - a)(y - b)\)
Exemplo:
Vamos aplicar a aproximação para a função \(f(x, y) = x^2 + y^2\) no ponto \((1, 1)\).
- Derivadas de Primeira Ordem:
\[ \dfrac{\partial f}{\partial x} = 2x \implies \dfrac{\partial f}{\partial x}(1, 1) = 2 \]
\[ \dfrac{\partial f}{\partial y} = 2y \implies \dfrac{\partial f}{\partial y}(1, 1) = 2 \]
- Derivadas de Segunda Ordem:
\[ \dfrac{\partial^2 f}{\partial x^2} = 2 \implies \dfrac{\partial^2 f}{\partial x^2}(1, 1) = 2 \]
\[ \dfrac{\partial^2 f}{\partial y^2} = 2 \implies \dfrac{\partial^2 f}{\partial y^2}(1, 1) = 2 \]
\[ \dfrac{\partial^2 f}{\partial x \partial y} = 0 \]
Substituindo na fórmula:
\[ f(x, y) \approx 2 + 2(x - 1) + 2(y - 1) + \dfrac{1}{2} \cdot 2 \cdot (x - 1)^2 + \dfrac{1}{2} \cdot 2 \cdot (y - 1)^2 \]
Simplificando:
\[ \begin{align} f(x, y) &\approx 2 + 2(x - 1) + 2(y - 1) + (x - 1)^2 + (y - 1)^2\\ f(x, y) &\approx (x - 1)^2 + (y - 1)^2 + 2x + 2y - 4 \end{align} \]
Vetor Gradiente
O vetor gradiente de uma função \(f(x, y)\) é um vetor composto pelas derivadas parciais de primeira ordem:
\[ \nabla f = \begin{bmatrix} \dfrac{\partial f}{\partial x} \quad\dfrac{\partial f}{\partial y} \end{bmatrix} \]
Ele aponta na direção de maior crescimento da função e sua magnitude indica a taxa de variação.
Matriz Hessiana
A matriz Hessiana contém as derivadas parciais de segunda ordem e tem a forma:
\[ \mathcal{H} = \begin{bmatrix} \dfrac{\partial^2 f}{\partial x^2} & \dfrac{\partial^2 f}{\partial x \partial y} \\ \dfrac{\partial^2 f}{\partial y \partial x} & \dfrac{\partial^2 f}{\partial y^2} \end{bmatrix} \]
Essa matriz descreve a curvatura da função.
Aproximação de Taylor de Segunda Ordem
Podemos escrever a fórmula da aproximação de Taylor de segunda ordem usando o vetor gradiente e a matriz Hessiana:
\[ f(x, y) \approx f(a, b) + \nabla f(a, b) \cdot \begin{bmatrix} x - a \\ y - b \end{bmatrix} + \dfrac{1}{2} \begin{bmatrix} x - a & y - b \end{bmatrix} \mathcal{H} \begin{bmatrix} x - a \\ y - b \end{bmatrix} \]
Decomposição da Fórmula:
Termo Constante: \(f(a, b)\)
Termo Linear (Gradiente):
\[ \nabla f(a, b) \cdot \begin{bmatrix} x - a \\ y - b \end{bmatrix} = \dfrac{\partial f}{\partial x}(a, b) \cdot (x - a) + \dfrac{\partial f}{\partial y}(a, b) \cdot (y - b) \]
- Termo Quadrático (Hessiana):
\[ \dfrac{1}{2} \begin{bmatrix} x - a & y - b \end{bmatrix} \mathcal{H} \begin{bmatrix} x - a \\ y - b \end{bmatrix} \]
Expandindo o produto matricial, temos:
\[ \dfrac{1}{2} \left[\dfrac{\partial^2 f}{\partial x^2}(x - a)^2 + 2 \dfrac{\partial^2 f}{\partial x \partial y}(x - a)(y - b) + \dfrac{\partial^2 f}{\partial y^2}(y - b)^2\right] \]
Exemplo:
Para a função \(f(x, y) = x^2 + y^2\) no ponto \((1, 1)\):
Gradiente:
\[ \nabla f = \begin{bmatrix} 2x \quad 2y \end{bmatrix} \implies \nabla f(1, 1) = \begin{bmatrix} 2 \quad 2 \end{bmatrix} \]
Hessiana:
\[ \mathcal{H} = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix} \]
Aproximação de Taylor de Segunda Ordem:
\[ f(x, y) \approx 2 + [2(x - 1) + 2(y - 1)] + \dfrac{1}{2} [2(x - 1)^2 + 2(y - 1)^2] \]
\[ \begin{align} f(x, y) &\approx 2 + 2(x - 1) + 2(y - 1) + (x - 1)^2 + (y - 1)^2\\ f(x, y) &\approx (x - 1)^2 + (y - 1)^2 + 2x + 2y - 4 \end{align} \]
A aproximação de Taylor de segunda ordem utiliza o vetor gradiente para representar a variação linear e a matriz Hessiana para descrever a curvatura da função. A combinação desses elementos permite uma aproximação mais precisa, levando em conta tanto a inclinação quanto a concavidade da superfície.
Determinante da hessiana
Conforme Goldman (2005), a curvatura de uma função com duas ou mais variáveis independentes, \(\kappa(x,y)\), é proporcional ao determinante da hessiana. O sinal da curvatura dependente apenas deste determinante.
\[ \kappa\propto \det(\mathcal{H}) \]
Exemplo: \(f(x,y)=e^{-x^2-y^2}\)
plot exp(-x^2-y^2) where x from -2 to 2 and y from -2 to 2
exp_neg_xy <- function(x, y) {
exp(-x^2 - y^2)
}
x_values <- seq(-2, 2, length.out = 100)
y_values <- seq(-2, 2, length.out = 100)
grid <- base::expand.grid(x = x_values, y = y_values)
grid$z <- base::with(grid, exp_neg_xy(x, y))
grid_melt <- reshape2::melt(grid, id.vars = c("x", "y"), measure.vars = "z")
ggplot_contour <- ggplot2::ggplot(grid_melt, ggplot2::aes(x = x, y = y, z = value)) +
ggplot2::geom_tile(ggplot2::aes(fill = value)) +
ggplot2::geom_contour(color = "black") +
ggplot2::scale_fill_gradientn(colours = grDevices::terrain.colors(10)) +
ggplot2::labs(title = "Gráfico de Contorno da Função e^{-x^2 - y^2}", x = "x", y = "y", fill = "z") +
ggplot2::theme_minimal()
print(ggplot_contour)
rayshader::plot_gg(ggplot_contour, multicore = TRUE, zoom = 0.8, phi = 40, theta = 30)
rayshader::render_snapshot(clear = FALSE)
fig <- plotly::plot_ly(x = ~x_values, y = ~y_values, z = ~matrix(grid$z, nrow = 100)) %>%
plotly::add_surface(
contours = list(
z = list(
show = TRUE,
usecolormap = TRUE,
highlightcolor = "#ff0000",
project = list(z = TRUE)
)
)
) %>%
plotly::layout(
scene = list(
xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "z"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = -0.64)
)
),
title = "Superfície 3D da Função e^{-x^2 - y^2}"
)
fig
f <- function(x, y) {
exp(-x^2 - y^2)
}
grad_f <- function(x, y) {
c(-2 * x * exp(-x^2 - y^2), -2 * y * exp(-x^2 - y^2))
}
x <- seq(-2, 2, length.out = 10)
y <- seq(-2, 2, length.out = 10)
grid <- expand.grid(x = x, y = y)
grid$z <- mapply(f, grid$x, grid$y)
gradients <- t(mapply(grad_f, grid$x, grid$y))
grid$u <- gradients[, 1]
grid$v <- gradients[, 2]
g <- ggplot2::ggplot(grid, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z), alpha = 0.3) +
ggplot2::geom_segment(ggplot2::aes(xend = x + u,
yend = y + v),
arrow = ggplot2::arrow(length = ggplot2::unit(0.3, "cm")),
color = "black") +
ggplot2::scale_fill_gradient(low = "white", high = "black") +
ggplot2::labs(title = "Vetor Gradiente de f(x, y) = exp(-x^2 - y^2)",
x = "x", y = "y") +
ggplot2::theme_minimal()
print(g)
A matriz Hessiana de \(e^{-x^2-y^2}\) é:
\[ \mathcal{H} = \exp(-(x^2 + y^2)) \begin{bmatrix} 4x^2 - 2 & 4xy \\ 4xy & 4y^2 - 2 \end{bmatrix} \]
O determinante da matriz Hessiana de \(e^{-x^2-y^2}\) é:
\[ \det(\mathcal{H})=4e^{-2(x^2 + y^2)}(1 - 2x^2 - 2y^2) \]
plot e^(-2 (x^2 + y^2)) (4 - 8 x^2 - 8 y^2) from x=-1 to x=1 from y=-1 to y=1 axes label "x" "y"
O sinal deste determinante depende apenas de \(1 - 2x^2 - 2y^2\).
A curva de inflexão, neste caso, é obtida por \(\det(\mathcal{H})=0\).
A solução é a equação da circunferência \(x^2 + y^2=\left(\frac{1}{\sqrt{2}}\right)^2\).
\(\det(\mathcal{H})>0\) implica \(x^2 + y^2<\left(\dfrac{1}{\sqrt{2}}\right)^2\).
x_values <- seq(-2, 2, length.out = 100)
y_values <- seq(-2, 2, length.out = 100)
grid <- expand.grid(x = x_values, y = y_values)
grid$condition <- (grid$x^2 + grid$y^2) < (1 / sqrt(2))^2
plot(grid$x, grid$y, col = ifelse(grid$condition, "darkgray", "white"),
pch = 15, cex = 0.6,
xlab = "x", ylab = "y",
main = expression(x^2 + y^2 < (frac(1, sqrt(2)))^2),
asp = 1)
theta <- seq(0, 2*pi, length.out = 200)
lines(1/sqrt(2) * cos(theta), 1/sqrt(2) * sin(theta), col = "black", lwd = 2)
abline(v=0,h=0,lty=2)
\(\det(\mathcal{H})<0\) implica \(x^2 + y^2>\left(\dfrac{1}{\sqrt{2}}\right)^2\).
x_values <- seq(-2, 2, length.out = 100)
y_values <- seq(-2, 2, length.out = 100)
grid <- expand.grid(x = x_values, y = y_values)
grid$condition <- (grid$x^2 + grid$y^2) > (1 / sqrt(2))^2
plot(grid$x, grid$y, col = ifelse(grid$condition, "darkgray", "white"),
pch = 15, cex = 0.6,
xlab = "x", ylab = "y",
main = expression(x^2 + y^2 > (frac(1, sqrt(2)))^2),
asp = 1)
theta <- seq(0, 2*pi, length.out = 200)
lines(1/sqrt(2) * cos(theta), 1/sqrt(2) * sin(theta), col = "black", lwd = 2)
abline(v=0,h=0,lty=2)
Exemplo: \(f(x,y)=-e^{-x^2-y^2}\)
plot -exp(-x^2-y^2) x from -2 to 2 and y from -2 to 2
exp_neg_xy <- function(x, y) {
-exp(-x^2 - y^2)
}
x_values <- seq(-2, 2, length.out = 100)
y_values <- seq(-2, 2, length.out = 100)
grid <- base::expand.grid(x = x_values, y = y_values)
grid$z <- base::with(grid, exp_neg_xy(x, y))
grid_melt <- reshape2::melt(grid, id.vars = c("x", "y"), measure.vars = "z")
ggplot_contour <- ggplot2::ggplot(grid_melt, ggplot2::aes(x = x, y = y, z = value)) +
ggplot2::geom_tile(ggplot2::aes(fill = value)) +
ggplot2::geom_contour(color = "black") +
ggplot2::scale_fill_gradientn(colours = grDevices::terrain.colors(10)) +
ggplot2::labs(title = "Gráfico de Contorno da Função e^{-x^2 - y^2}", x = "x", y = "y", fill = "z") +
ggplot2::theme_minimal()
print(ggplot_contour)
rayshader::plot_gg(ggplot_contour, multicore = TRUE, zoom = 0.8, phi = 40, theta = 30)
rayshader::render_snapshot(clear = FALSE)
fig <- plotly::plot_ly(
x = ~x_values,
y = ~y_values,
z = ~matrix(grid$z, nrow = 100)
)
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 = "z"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = -0.64)
)
),
title = "Superfície 3D da Função e^{-x^2 - y^2}"
)
fig
# Definir a nova função
f <- function(x, y) {
-exp(-x^2 - y^2)
}
# Definir o gradiente da função
grad_f <- function(x, y) {
c(2 * x * exp(-x^2 - y^2), 2 * y * exp(-x^2 - y^2))
}
# Gerar uma grade de pontos
x <- seq(-2, 2, length.out = 10)
y <- seq(-2, 2, length.out = 10)
grid <- expand.grid(x = x, y = y)
# Calcular os valores da função e do gradiente
grid$z <- mapply(f, grid$x, grid$y)
gradients <- t(mapply(grad_f, grid$x, grid$y))
grid$u <- gradients[, 1]
grid$v <- gradients[, 2]
# Criar o gráfico
library(ggplot2)
g <- ggplot2::ggplot(grid, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z), alpha = 0.3) +
ggplot2::geom_segment(ggplot2::aes(xend = x + u,
yend = y + v),
arrow = ggplot2::arrow(length = ggplot2::unit(0.3, "cm")),
color = "black") +
ggplot2::scale_fill_gradient(low = "white", high = "black") +
ggplot2::labs(title = "Vetor Gradiente de f(x, y) = -exp(-x^2 - y^2)",
x = "x", y = "y") +
ggplot2::theme_minimal()
print(g)
A matriz Hessiana de \(-e^{-x^2-y^2}\) é:
\[ \mathcal{H} = \exp(-(x^2 + y^2)) \begin{bmatrix} 2-4x^2 & -4xy \\ -4xy & 2-4y^2 \end{bmatrix} \]
O determinante da matriz Hessiana de \(e^{-x^2-y^2}\) é:
\[ \det(\mathcal{H})=4e^{-2(x^2 + y^2)}(1 - 2x^2 - 2y^2) \]
O determinante da hessiana das duas funções é o mesmo. Portanto, o ponto extremante (mínimo ou máximo) local pode ser encontrado na região com determinante da hessiana positivo, \(\det(\mathcal{H})>0\).
second derivative exp(-(x^2+y^2)) for x
Hessian Exp(-x^2-y^2) = Hessian -Exp(-x^2-y^2)
[1] "-(exp(-x^2 - y^2) * (2 * x))" "-(exp(-x^2 - y^2) * (2 * y))" [,1]
[1,] "-(exp(-x^2 - y^2) * 2 - exp(-x^2 - y^2) * (2 * x) * (2 * x))"
[2,] "exp(-x^2 - y^2) * (2 * x) * (2 * y)"
[,2]
[1,] "exp(-x^2 - y^2) * (2 * y) * (2 * x)"
[2,] "-(exp(-x^2 - y^2) * 2 - exp(-x^2 - y^2) * (2 * y) * (2 * y))"[1] 0 0 [,1] [,2]
[1,] -2 0
[2,] 0 -2[1] 4[1] -0.2706706 -0.2706706 [,1] [,2]
[1,] 0.2706706 0.5413411
[2,] 0.5413411 0.2706706[1] -0.2197877[1] -0.019996 -0.019996 [,1] [,2]
[1,] -1.99920012 0.00039992
[2,] 0.00039992 -1.99920012[1] 3.996801[1] "exp(-x^2 - y^2) * (2 * x)" "exp(-x^2 - y^2) * (2 * y)" [,1]
[1,] "exp(-x^2 - y^2) * 2 - exp(-x^2 - y^2) * (2 * x) * (2 * x)"
[2,] "-(exp(-x^2 - y^2) * (2 * x) * (2 * y))"
[,2]
[1,] "-(exp(-x^2 - y^2) * (2 * y) * (2 * x))"
[2,] "exp(-x^2 - y^2) * 2 - exp(-x^2 - y^2) * (2 * y) * (2 * y)"[1] 0 0 [,1] [,2]
[1,] 2 0
[2,] 0 2[1] 4[1] 0.2706706 0.2706706 [,1] [,2]
[1,] -0.2706706 -0.5413411
[2,] -0.5413411 -0.2706706[1] -0.2197877[1] 0.019996 0.019996 [,1] [,2]
[1,] 1.99920012 -0.00039992
[2,] -0.00039992 1.99920012[1] 3.996801Otimização de função de duas variáveis (sem restrições)
Ponto crítico e extremante
Extremante é um ponto de máximo ou de mínimo local da função.
Ponto crítico é um ponto que anula o gradiente.
Extremante pode ser ponto crítico, mas nem todo ponto crítico é extremante.
Fig. As derivadas parciais são nulas no ponto extremante local.
Um ponto de sela é um ponto crítico de uma função onde a função não possui nem um máximo nem um mínimo local, mas sim uma mudança de concavidade. O nome vem da semelhança com a forma de uma sela de cavalo: curva para cima em uma direção e para baixo em outra.
Fig. Ponto de sela de z = y22
Condição necessária
A condição necessária para extremante é a existência de ponto crítico, i.e., \(\nabla f\left(x^*, y^*\right)=0\).
Fig. 12.3. Encontrando o máximo de função.
\(f(x, y) = \sqrt{xy}\)
O gradiente de \(f\) é dado por:
\[ \nabla f(x, y) = \left( \dfrac{\partial f}{\partial x}, \dfrac{\partial f}{\partial y} \right) \]
Calculando as derivadas parciais:
\[ \dfrac{\partial f}{\partial x} = \dfrac{y}{2\sqrt{xy}} \]
\[ \dfrac{\partial f}{\partial y} = \dfrac{x}{2\sqrt{xy}} \]
Portanto, o gradiente de \(f(x, y) = \sqrt{xy}\) é:
\[ \nabla f(x, y) = \left( \dfrac{y}{2\sqrt{xy}}, \dfrac{x}{2\sqrt{xy}} \right) \]
Note que o gradiente não pode ser nulo. Portanto, não há ponto crítico e, desta forma, não há extremante.
\(g(x, y) = \exp(-x^2 - y^2)\)
O gradiente de \(g\) é dado por:
\[ \nabla g(x, y) = \left( \dfrac{\partial g}{\partial x}, \dfrac{\partial g}{\partial y} \right) \]
Calculando as derivadas parciais:
\[ \dfrac{\partial g}{\partial x} = -2x \exp(-x^2 - y^2) \]
\[ \dfrac{\partial g}{\partial y} = -2y \exp(-x^2 - y^2) \]
Portanto, o gradiente é:
\[ \nabla g(x, y) = \exp(-x^2 - y^2)\left( -2x , -2y \right) \]
\[ \nabla g(0, 0) =\exp(-0^2 - 0^2)\left( -2\times 0 , -2\times 0 \right)=(0,0) \]
Note que o ponto crítico é \((0,0)\). Ainda não temos condição de determinar se este ponto crítico é extremante.
Condição suficiente
A condição suficiente de segunda ordem de ponto extremante demonstrada em Teichroew (1964) usando principal minor e está relacionada com a teoria de autovalor (eigenvalue) de álgebra linear (Moylan, 1977; Venit & Katz, 2000; González-Vega et al., 2023).
A condição suficiente de segunda ordem determina se um ponto crítico é extremante, i.e., de mínimo ou máximo local é:
\[ \begin{align} \det(\mathcal{H}(f\left(x^*, y^*\right)))&=f_{xx}\left(x^*, y^*\right)f_{yy}\left(x^*, y^*\right) - \left(f_{xy}\left(x^*, y^*\right)\right)^2\\ \det(\mathcal{H}(f\left(x^*, y^*\right)))&>0 \end{align} \]
O ponto crítico \(\left(a^*, b^*\right)\) da função \(f(x,y)\) é ponto de mínimo local se:
\[ \begin{align} \det(\mathcal{H}(f\left(x^*, y^*\right)))&>0\\ f_{xx}\left(x^*, y^*\right)&>0 \end{align} \]
Conforme Neuhauser & Roper, 2018, p. 614, a matriz Hessiana é definida positiva: no caso geral, há autovalores apenas positivos; neste caso, os produtórios acumulados dos autovalores têm sinais positivos.
O ponto crítico \(\left(a^*, b^*\right)\) da função \(f(x,y)\) é ponto de máximo local se:
\[ \begin{align} \det(\mathcal{H}(f\left(x^*, y^*\right)))>0\\ f_{xx}\left(x^*, y^*\right)<0 \end{align} \]
Conforme Neuhauser & Roper, 2018, p. 614, matriz Hessiana é definida negativa: no caso geral, há autovalores apenas negativos; neste caso, os produtórios acumulados dos autovalores têm sinais alternados.
O ponto crítico \(\left(a^*, b^*\right)\) da função \(f(x,y)\) não é ponto de mínimo ou máximo, se:
\[ \det(\mathcal{H}(f\left(x^*, y^*\right)))<0 \]
Conforme Neuhauser & Roper, 2018, p. 614, a matriz Hessiana é indefinida: no caso geral, há autovalores positivo(s) e negativo(s); neste caso, os produtórios acumulados dos autovalores têm sinais não definidos.
O ponto crítico \(\left(a^*, b^*\right)\) da função \(f(x,y)\) pode ser ponto de mínimo ou máximo se:
\[ \det(\mathcal{H}(f\left(x^*, y^*\right)))=0 \]
Conforme Neuhauser & Roper, 2018, p. 614, a matriz Hessiana é semi-definida: no caso geral, há autovalores positivo(s) e nulo(s); o produtório de todos os autovalores é nulo; note que zero não é um número positivo.
\(f(x, y) = \sqrt{xy}\)
Gradiente:
\[ \nabla f(x, y) = \left( \dfrac{y}{2\sqrt{xy}}, \dfrac{x}{2\sqrt{xy}} \right) \]
Não há ponto crítico.
Hessiana:
\[ \mathcal{H}_f = \begin{pmatrix} -\dfrac{y^2}{4 (x y)^{3/2}} & \dfrac{xy}{4 (x y)^{3/2}} \\ \dfrac{xy}{4 (x y)^{3/2}} & -\dfrac{x^2}{4 (x y)^{3/2}} \end{pmatrix} \]
\(g(x, y) = \exp(-x^2 - y^2)\)
Gradiente:
\[ \nabla g(x, y) = \exp(-x^2 - y^2)\left[ -2x , -2y \right] \]
Ponto crítico: \((0,0)\)
Hessiana:
\[ \mathcal{H}_g = \exp(-x^2 - y^2) \begin{bmatrix} 4x^2 - 2 & 4xy \\ 4xy & 4y^2 - 2 \end{bmatrix} \]
Determinante da Hessiana:
\[ \det(\mathcal{H}_g) = -8 \exp\left(-2\left(x^2 + y^2\right)\right) \left(x^2 + y^2 - \frac{1}{2}\right) \]
Hessiana em \((0,0)\):
\[ \mathcal{H}_g = \begin{bmatrix} -2 & 0 \\ 0 & -2 \end{bmatrix} \]
Determinante da Hessiana em \((0,0)\):
\[ \det(\mathcal{H}_g) = 4 \]
O ponto crítico é extremante de máximo.
\(h(x, y) = -\exp(-x^2 - y^2)\)
Gradiente:
\[ \nabla h(x, y) = \exp(-x^2 - y^2)\left[2x , 2y \right] \]
Ponto crítico: \((0,0)\)
Hessiana:
\[ \mathcal{H}_h = \exp(-x^2 - y^2) \begin{bmatrix} -(4x^2 - 2) & -4xy \\ -4xy & -(4y^2 - 2) \end{bmatrix} \]
Determinante da Hessiana:
\[ \det(\mathcal{H}_g) = -8 \exp\left(-2\left(x^2 + y^2\right)\right) \left(x^2 + y^2 - \frac{1}{2}\right) \]
Hessiana em \((0,0)\):
\[ \mathcal{H}_h = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix} \]
Determinante da Hessiana em \((0,0)\):
\[ \det(\mathcal{H}_g) = 4 \]
O ponto crítico é extremante de mínimo.
Exemplo: Transporte fluvial de produto granular
Quatrocentos metros cúbicos de um produto granular devem ser transportados para outro lado do rio por meio de um serviço de ferryboat. Qualquer quantidade do produto pode ser transportada em uma viagem que custa dez centavos por viagem. Os ancoradouros do ferryboat estão construídos de modo que o contêiner no qual o produto vai ser transportado deve ter meio metro de profundidade, ao passo que a largura e o comprimento podem variar. Os custos por metro quadrado das superfícies terminais (da parte da frente à parte detrás), laterais e do fundo do contêiner são R$ 20, R$ 10 e R$ 10, respectivamente. Se o tamanho do contêiner for pequeno, seu custo poderá ser baixo, mas o número de viagens será maior, aumentando o custo do transporte e vice-versa. Quais são os valores do comprimento e largura em metro do contêiner, isto é, quais as dimensões do contêiner que minimizam o custo total do transporte do produto granular?
Dimensões do contêiner.
Solução:
- \(a\): largura da base retangular do contêiner (m)
- \(b\): comprimento da base retangular do contêiner (m)
- Custo das superfícies terminais: \(20a\dfrac{1}{2}2=20a\) (R$)
- Custo das superfícies laterais: \(10b\dfrac{1}{2}2=10b\) (R$)
- Custo da superfície do fundo: \(10ab\) (R$)
- Quantidade de viagens: \(\dfrac{400}{\dfrac{1}{2}ab}=\dfrac{800}{ab}\) (viagem)
- Custo da quantidade de viagens: \(0.10\dfrac{400}{\dfrac{1}{2}ab}=\dfrac{80}{ab}\) (R$)
- \(C(a,b)\): custo total do transporte do produto granular (R$).
A função do custo total de transporte do produto granular por via fluvial é dada por:
\[ C(a,b)=\dfrac{80}{ab}+20a+10b+10ab \]
O contêiner de 0.5 m de altura, 1 m de largura e 2 m de comprimento, isto é, de 1 m3, realiza 400 viagens. O custo total mínimo do transporte fluvial de produto granular é R$ 100, isto é, \(C_{\text{min}}(1,2)=100\).
f <- function(a, b) {
80 / (a * b) + 20 * a + 10 * b + 10 * a * b
}
a_values <- seq(0.5, 2, length.out = 100)
b_values <- seq(1, 3, length.out = 100)
grid <- base::expand.grid(a = a_values, b = b_values)
grid$z <- base::with(grid, f(a, b))
grid_melt <- reshape2::melt(grid,
id.vars = c("a", "b"),
measure.vars = "z")
ggplot_contour <- ggplot2::ggplot(grid_melt,
ggplot2::aes(x = a,
y = b,
z = value)) +
ggplot2::geom_tile(ggplot2::aes(fill = value)) +
ggplot2::geom_contour(color = "black") +
ggplot2::scale_fill_gradientn(colours = grDevices::terrain.colors(10)) +
ggplot2::labs(title = "Gráfico de Contorno",
x = "a", y = "b", fill = "C") +
ggplot2::theme_minimal()
print(ggplot_contour)
rayshader::plot_gg(ggplot_contour,
multicore = TRUE,
zoom = 0.8, phi = 40, theta = 30)
rayshader::render_snapshot(clear = FALSE)
fig <- plotly::plot_ly(x = ~a_values,
y = ~b_values,
z = ~base::matrix(grid$z, nrow = 100))
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 = "a"),
yaxis = list(title = "b"),
zaxis = list(title = "C"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = -0.64)
)
),
title = "Superfície 3D da Função")
figgrad_f <- function(a, b) {
c(-80 / (a^2 * b) + 20 + 10 * b, -80 / (a * b^2) + 10 + 10 * a)
}
a <- seq(0.5, 2, length.out = 5)
b <- seq(1, 3, length.out = 3)
grid <- expand.grid(a = a, b = b)
grid$z <- mapply(f, grid$a, grid$b)
gradients <- t(mapply(grad_f, grid$a, grid$b))
grid$u <- gradients[, 1]
grid$v <- gradients[, 2]
g <- ggplot2::ggplot(grid, ggplot2::aes(x = a, y = b)) +
ggplot2::geom_tile(ggplot2::aes(fill = z), alpha = 0.3) +
ggplot2::geom_segment(ggplot2::aes(xend = a + u,
yend = b + v),
arrow = ggplot2::arrow(length =
ggplot2::unit(0.2,
"cm")),
color = "darkgray") +
ggplot2::scale_fill_gradient(low = "white", high = "black") +
ggplot2::labs(title = "Vetor Gradiente de f(a, b) = 80/(ab) + 20a + 10b + 10ab",
x = "a", y = "b") +
ggplot2::geom_point(ggplot2::aes(x = 1, y = 2), color = "blue",
size = 3) +
ggplot2::theme_minimal()
print(g)
extrema 80/(a b) + 20 a+ 10 b + 10 a b for a>0 and b>0
minimize 80/(a b) + 20 a+ 10 b + 10 a b for a>0 and b>0
Plot[80/ab+20a+10b+10ab, {a,0.5,2}, {b,1,3}]
A solução exata do problema do custo total mínimo do transporte fluvial de produto granular é dada pelas condições de primeira e segunda ordens. A condição necessária de primeira ordem consiste na determinação dos pontos críticos, isto é, candidatos a ponto extremante de mínimo, por meio da resolução do sistema de equações de duas derivadas parciais da função de custo total, ou seja:
\[ C_{a}\left(a^*, b^*\right)=0 \\ C_{b}\left(a^*, b^*\right)=0 \]
A primeira equação \(C_{a}\left(a^*, b^*\right)=0\) contém no lado esquerdo uma derivada parcial da função do custo total em relação à dimensão \(a\) do contêiner.
A segunda equação \(C_{b}\left(a^*, b^*\right)=0\) contém no lado esquerdo uma derivada parcial da função do custo total em relação à dimensão \(b\) do contêiner.
A ideia subjacente ao sistema de derivadas parciais igualadas a zero é que há uma superfície plana paralela ao plano do domínio que tangencia a superfície convexa do custo total no ponto de mínimo.
Portanto:
\[ C_{a}\left(a^*, b^*\right)=-\dfrac{80}{(a^{*})^2b^*}+20+10b^*=0 \\ C_{b}\left(a^*, b^*\right)=-\dfrac{80}{a^{*}(b^{*})^2}+10+10a^{*}=0 \]
A solução do sistema de equações é \(a^*=1\) m e \(b^*=2\) m.
(D[80/(a b) + 20 a +10 b + 10 a b, a], D[80/(a b)+20 a+10 b+10 a b, b]) = (0, 0)
Portanto, há um ponto crítico no problema, isto é, \(\left(a^*, b^*\right)=(1,2)\).
Resta saber se o ponto crítico é de mínimo ou máximo da função.
Um ponto crítico pode ser de mínimo (máximo) se neste ponto a função é convexa (côncava). O tipo de concavidade da superfície da função no ponto depende de sua curvatura (Goldman, 2005).
Se a função é convexa (côncava) na região e tem extremante, ele é um ponto de mínimo (máximo).
O extremante \(\left(x^*, y^*\right)\) é mínimo, se:
\[ \begin{align} f_{xx}\left(x^*, y^*\right)f_{yy}\left(x^*, y^*\right) - \left(f_{xy}\left(x^*, y^*\right)\right)^2>0\\ f_{xx}\left(x^*, y^*\right)>0 \end{align} \]
f <- expression(80/(a*b) + 20*a + 10*b + 10*a*b)
gradiente <- calculus::gradient(f, var=c("a","b"))
gradiente[1] "20 - 80 * b/(a * b)^2 + 10 * b" "10 - 80 * a/(a * b)^2 + 10 * a" [,1]
[1,] "80 * b * (2 * (b * (a * b)))/((a * b)^2)^2"
[2,] "10 - (80/(a * b)^2 - 80 * a * (2 * (b * (a * b)))/((a * b)^2)^2)"
[,2]
[1,] "10 - (80/(a * b)^2 - 80 * b * (2 * (a * (a * b)))/((a * b)^2)^2)"
[2,] "80 * a * (2 * (a * (a * b)))/((a * b)^2)^2" [1] 0 0 [,1] [,2]
[1,] 80 30
[2,] 30 20[1] 700No problema, tem-se:
\[ \begin{align} C_{aa}\left(a^*, b^*\right)C_{bb}\left(a^*, b^*\right)-\left(C_{ab}\left(a^*, b^*\right)\right)^2&=\dfrac{25600}{(a^*b^*)^4}-\left(\dfrac{80}{(a^*b^*)^2}+10\right)^2=700\\ C_{aa}\left(a^*, b^*\right)&=\dfrac{160}{(a^*)^3b^*}=80 \end{align} \]
Portanto, o ponto crítico \(\left(a^*, b^*\right)=(1,2)\) é de mínimo local. Como a função é convexa em todo domínio, então o ponto é de mínimo global.
D[80/(a b)+20 a + 10 b + 10 a b, a, a]
D[80/(a b)+20 a + 10 b + 10 a b, b, b]
D[80/(a b)+20 a + 10 b + 10 a b, a, b]
Generalizando a otimização
Para calcular os principal minors da matriz Hessiana fornecida, precisamos calcular os determinantes das submatrizes principais. A matriz Hessiana é uma matriz simétrica quadrada, e os principal minors são os determinantes das submatrizes obtidas ao remover as mesmas colunas e linhas da matriz original.
Dada a matriz Hessiana:
\[ \mathcal{H} = \begin{pmatrix} 80 & 30 \\ 30 & 20 \end{pmatrix} \]
Vamos calcular os determinantes das submatrizes principais:
- Determinante da matriz \(1 \times 1\) obtida ao considerar apenas o primeiro elemento:
\[ \mathcal{H}_1 = \begin{pmatrix} 80 \end{pmatrix} \implies \det(\mathcal{H}_1) = 80 \]
- Determinante da matriz \(2 \times 2\) original:
\[ \mathcal{H}_2 = \begin{pmatrix} 80 & 30 \\ 30 & 20 \end{pmatrix} \implies \det(\mathcal{H}_2) = 80 \cdot 20 - 30 \cdot 30 = 1600 - 900 = 700 \]
Portanto, os principal minors da matriz Hessiana são 80 e 700.
Para determinar a natureza dos pontos críticos usando a matriz Hessiana, devemos verificar os sinais dos seus principal minors. Se todos os principal minors tiverem sinais positivos e não alternantes, isso indica que o ponto crítico é um mínimo local.
No nosso caso, calculamos os principal minors da matriz Hessiana:
\[ \mathcal{H}=\begin{pmatrix} 80 & 30 \\ 30 & 20 \end{pmatrix} \]
Ambos os principal minors são positivos:
- \(\det(\mathcal{H}_1) = 80 > 0\)
- \(\det(\mathcal{H}_2) = 700 > 0\)
Como ambos são positivos e não alternam de sinal, isso indica que a matriz Hessiana é definida positiva. Portanto, o ponto crítico associado a essa matriz Hessiana é um mínimo local extremo.
Matriz definida positiva tem autovalores positivos. Neste caso, a matriz Hessiana tem autovalores aproximadamente iguais a 92.43 e 7.57.
A análise correta dos sinais dos principal minors para determinar a natureza de um ponto crítico é dada pelos seguintes critérios:
- Mínimo Local: Os dois principal minors da matriz Hessiana são positivos: \(\det(\mathcal{H}_1)>0\) e \(\det(\mathcal{H}_2)>0\).
- Máximo Local: Os dois principal minors alternam de sinal começando com um sinal negativo: \(\det(\mathcal{H}_1)<0\) e \(\det(\mathcal{H}_2)>0\).
- Ponto de Sela: O segundo principal minor é negativo: \(\det(\mathcal{H}_2)<0\).
Dado os principal minors calculados:
\[ \begin{align} \mathcal{H}_1 &= 80 \quad (\text{positivo})\\ \mathcal{H}_2 &= 700 \quad (\text{positivo}) \end{align} \]
Ambos são positivos, o que indica que a matriz Hessiana é definida positiva. Portanto, o ponto crítico é um mínimo local.
# Definir a matriz Hessiana
hessian_matrix <- matrix(c(80, 30, 30, 20), nrow = 2, byrow = TRUE)
# Calcular os autovalores
eigenvalues <- eigen(hessian_matrix)$values
# Calcular os principal minors
principal_minor_1 <- det(hessian_matrix[1, 1, drop = FALSE])
principal_minor_2 <- det(hessian_matrix)
# Determinar a natureza do ponto crítico
if (principal_minor_1 > 0 && principal_minor_2 > 0) {
point_nature <- "Mínimo Local"
}
if (principal_minor_1 < 0 && principal_minor_2 > 0) {
point_nature <- "Máximo Local"
}
if (principal_minor_2 < 0) {
point_nature <- "Ponto de Sela"
}
if (principal_minor_2 == 0) {
point_nature <- "Indeterminado"
}
# Exibir os resultados
cat("Autovalores:\n")Autovalores:[1] 92.426407 7.573593Principal minors:Principal minor 1: 80 Principal minor 2: 700 Natureza do ponto crítico: Mínimo Local A seguir, o critério para \(n\ge3\) variáveis independentes é generalizado.
Seja \((x_1^{\ast}, x_2^{\ast}, \ldots,x_n^{\ast})\) um ponto crítico.
Este ponto crítico é ponto de mínimo local se todos os autovalores da matriz Hessiana são positivos ou, equivalentemente:
\[ \det(\mathcal{H}_1)>0, \det(\mathcal{H}_2)>0, \ldots,\det(\mathcal{H}_n)>0 \]
Este ponto crítico é ponto de máximo local se todos os autovalores da matriz Hessiana são negativos ou, equivalentemente:
\[\det(\mathcal{H}_1)<0, \det(\mathcal{H}_2)>0, \det(\mathcal{H}_3)<0,\ldots\]
simplify hessian matrix x^2+y^2+z^2
simplify hessian matrix -(x^2+y^2+z^2)
simplify hessian matrix x^2+y^2-z^2
simplify hessian matrix x^2+y^2-z^2
saddle point - x^2 - y^2 - z^2 + y z + x + 2 z
extrema - x^2 - y^2 - z^2 + y z + x + 2 z
simplify hessian matrix - x^2 - y^2 - z^2 + y z + x + 2 z
Exemplo: Ajustamento de curva pelo método de mínimos quadrados
Há fenômenos nas ciências biológicas que não permitem a adoção de uma função formulada matematicamente a partir de suposições teóricas. Às vezes, é necessário estimar a função a partir de dados observados. O fenômeno mais simples é aquele no qual os dados observados consistem em uma coleção finita de pares ordenados, e uma coordenada é o valor da variável dependente e a outra é o valor da variável independente.
A coleção de pontos é ordenada pelos valores da variável independente e o gráfico pode ser desenhado, sendo essa variável uma abscissa. Uma condição importante é que, para cada valor observado da variável independente, há apenas um valor observado da variável dependente. Surge então a pergunta: qual é a função que melhor representa os pontos observados no domínio observado. O domínio observado consiste no intervalo entre o menor e o maior valores observados da variável independente.
Há duas formas de responder a essa questão. A primeira abordagem consiste em determinar os parâmetros de uma única função que melhor representa a disposição dos pontos no gráfico sem necessariamente passar por todos os pontos observados. Essa abordagem é denominada ajustamento de curva. A segunda abordagem consiste em ligar os pontos por meio de diferentes polinômios de terceiro grau, isto é, por splines cúbicos ou naturais. Essa abordagem é denominada interpolação natural ou por spline cúbico.
O intuito das duas abordagens é preencher os valores ausentes da variável dependente dentro do domínio observado. Evidentemente, se há apenas dois pontos distintos observados, o melhor modelo de função é o afim. A função adotada para ajustamento de curva possui necessariamente parâmetros que precisam ser estimados a partir dos dados observados. Uma competição entre funções plausíveis para modelar o fenômeno pode ocorrer e, dessa forma, é necessário estabelecer um critério para selecionar a melhor função estimada.
Além das funções polinomiais, as curvas de crescimento também podem ser utilizadas para modelar. Os valores da variável independente não precisam ser equiespaçados. Do ponto de vista artístico, as curvas francesas podem ser usadas para desenhar uma curva contínua e não lisa para interligar os pontos. A interpolação natural é a que produz a curva que menos dista da curva original se os pontos são alterados de posição, isto é, ela é a mais robusta diante de modificações nos dados. Além disso, a curva gerada por interpolação natural ou por splines cúbicos é a que minimiza a curvatura total no domínio observado, ou seja, a curva com a tensão máxima que liga os pontos mantendo a continuidade (sem salto) e suavidade (sem cúspide) – um formato orgânico da função.
A interpolação natural resulta em curvas cúbicas ligando os pontos observados. O ajustamento de curva difere da interpolação pelo fato de exigir apenas uma função polinomial que melhor represente os pontos sem necessariamente passar por eles. Na interpolação exige-se que as funções resultantes passem pelos pontos observados. No ajustamento de curva não se exige que a função resultante passe pelos pontos observados. Isso pode ocorrer, mas não é compulsório.
Dessa maneira, no ajustamento de curva é necessário experimentar algumas formas de funções e escolher a que melhor represente os pontos observados. Naturalmente, se houver apenas dois pontos distintos observados, então a curva de melhor ajustamento é a retilínea.
Ajustar uma curva linear, quadrática e cúbica e selecionar aquela com melhor ajustamento aos dados observados da tabela a seguir:
| \(i\) | \(t_i\) (ano) | \(Y_i\) (milhão de R$/ano) |
|---|---|---|
| 1 | 1 | 5 |
| 2 | 2 | 8 |
| 3 | 3 | 25 |
Modelo linear
\[ Y_i=\beta_0+\beta_1t_i+\varepsilon_i\\ i=1,2,3 \]
O modelo linear tem os seguintes elementos:
- variável dependente \(Y\),
- variável independente \(t\),
- parâmetros desconhecidos \(\beta_0\) e \(\beta_1\)
- termo de erro \(\varepsilon\).
O domínio observado é \(t\in[1,3]\).
O termo de erro \(\varepsilon\) é o valor que adicionado à curva retilínea resulta no valor observado da variável dependente \(Y\).
O sistema linear nos parâmetros é o seguinte:
\[ \begin{align} 5&=\beta_0+\beta_1+\varepsilon_1\\ 8&=\beta_0+2\beta_1+\varepsilon_2\\ 25&=\beta_0+3\beta_1+\varepsilon_3 \end{align} \]
Observe que há 2 parâmetros e 3 termos de erro desconhecidos totalizando 5 incógnitas a serem estimadas com 3 pontos.
Seriam necessárias mais 2 equações não redundantes para estimar de modo único as 5 incógnitas. No entanto, essas duas equações adicionais não existem.
É preciso impor o critério dos mínimos quadrados ordinários (OLS), i.e., de minimização da soma dos quadrados dos termos de erro.
Matematicamente, a quantidade a ser minimizada em relação aos dois parâmetros pelos mínimos quadrados é dada pela soma de quadrados dos termos de erro:
\[ Q=\sum_{i=1}^{3}{\varepsilon_{i}^{2}} \]
As variáveis independentes da função anterior \(Q\) são os parâmetros \(\beta_0\) e \(\beta_1\).
O lado direito da equação deve ser expresso em função dos parâmetros \(\beta_0\) e \(\beta_1\) para a aplicação das condições de primeira e de segunda ordens.
Os termos de erro são expressos em função dos parâmetros a serem estimados:
\[ \varepsilon_i=Y_i-\beta_0-\beta_1t_i\\ i=1,2,3 \]
Dessa forma, tem-se que:
\[ Q=\sum_{i=1}^{3}\left({Y_i-\beta_0-\beta_1t_i}\right)^2 \]
A condição de primeira ordem para a determinação do ponto crítico é dada pelas derivadas parciais da soma dos quadrados dos termos de erro em relação a cada um dos parâmetros igualadas a zero:
\[ \begin{align} Q_{\beta_0}&=-2\sum_{i=1}^{3}\left({Y_i-b_0-b_1t_i}\right)=0\\ Q_{\beta_1}&=-2\sum_{i=1}^{3}\left({Y_i-b_0-b_1t_i}\right)t_i=0 \end{align} \]
O ponto crítico \((b_0,b_1)\) é obtido pela resolução do sistema de equações de derivadas parciais anterior e seus estimadores pontuais são:
\[ \begin{align} b_1&=\dfrac{\sum_{i=1}^{3}{\left(t_i-\dfrac{\sum_{i=1}^{3}{t_i}}{3}\right)Y_i}}{\sum_{i=1}^{3}{\left(t_i-\dfrac{\sum_{i=1}^{3}{t_i}}{3}\right)^2}}=\dfrac{\sum\left(t_i-\bar{t}\right)Y_i}{\sum{\left(t_i-\bar{t}\right)^2}}\\ b_0&=\dfrac{\sum_{i=1}^{3}{Y_i}}{3}-b_1\dfrac{\sum_{i=1}^{3}{t_i}}{3}=\bar{Y}-b_1\bar{t} \end{align} \]
Os valores estimados dos parâmetros são os seguintes:
\[ \begin{align} b_1&=\dfrac{(1-2)5+(2-2)8+(3-2)25}{(1-2)^2+(2-2)^2+(3-2)^2}=10\\ b_0&=\dfrac{5+8+25}{3}-10\dfrac{1+2+3}{3}=-7.\bar{3} \end{align} \]
O extremante \((b_0,b_1)=(-7.\bar{3},10)\) minimiza a soma de quadrados dos termos de erro, pois as seguintes desigualdades são sempre válidas:
\[ \begin{align} Q_{\beta_0\beta_0}\left(b_0,b_1\right)Q_{\beta_1\beta_1}\left(b_0,b_1\right)-\left(Q_{\beta_0\beta_1}\left(b_0,b_1\right)\right)^2>0\\ Q_{\beta_0\beta_0}\left(b_0,b_1\right)>0\\ \end{align} \]
No problema, temos:
\[ \begin{align} Q_{\beta_0\beta_0}\left(b_0,b_1\right)Q_{\beta_1\beta_1}\left(b_0,b_1\right)-\left(Q_{\beta_0\beta_1}\left(b_0,b_1\right)\right)^2&=\left(\sum_{i=1}^{3}{2}\right) \left(2\sum_{i=1}^{3}{t_{i}^{2}}\right)-\left(2\sum_{i=1}^{3}{t_{i}}\right)^2=12\\ Q_{\beta_0\beta_0}\left(b_0,b_1\right)&=\sum_{i=1}^{3}{2}=6 \end{align} \]
Portanto, o modelo linear estimado por OLS é:
\[ y=-7.\bar{3}+10t\\ t\in[1,3] \]
Note que \(t\) e \(y\) são variáveis contínuas (reais) no domínio observado.
O resíduo é a estimativa do termo de erro:
\[ e_i=Y_i-y_i\\ i=1,2,3 \]
A soma dos dos resíduos é sempre nula:
\[ \sum_{i=1}^{3}{e_{i}}=0 \]
No modelo linear, a reta passa sempre pelo centróide:
\[ y\left(\bar{t}\right)=\bar{Y} \]
LeastSquares[{{1, 1}, {1, 2}, {1, 3}}, {5, 8, 25}]
linear fit {{1,5},{2,8},{3,25}}
Fit[{{1, 5}, {2, 8}, {3, 25}}, {1, t}, t]
FindFit[{{1, 5}, {2, 8}, {3, 25}}, \beta_0+\beta_1 t, {\beta_0, \beta_1}, t]
NonlinearModelFit[{{1, 5}, {2, 8}, {3, 25}}, \beta_0+\beta_1 t, {\beta_0, \beta_1}, t]
A notação matricial permite resolver o problema e possibilita ganho de escala para ajustamento de curva polinomial de ordem maior que 1.
A fórmula geral do método OLS é dada por:
\[ \Large b=\left(x^Tx\right)^{-1}x^TY \]
Sendo que \(x^T\) é a matriz transposta da matriz de planejamento (design matrix: Wikipedia) \(x\).
A matriz \(x^Tx\) é conhecida como matriz de informação. A matriz de informação é quadrada e simétrica.
\[ \begin{align} b&=\left[\begin{array}{r} b_0\\ b_1 \end{array}\right] = \left[\begin{array}{r} -7.\bar{3}\\ 10 \end{array}\right]\\ Y&=\left[\begin{array}{r} Y_1\\ Y_2\\ Y_3 \end{array}\right] = \left[\begin{array}{r} 5\\ 8\\ 25 \end{array}\right]\\ x&=\left[\begin{array}{rr} 1 & t_1\\ 1 & t_2\\ 1 & t_3 \end{array}\right] = \left[\begin{array}{rr} 1 & 1\\ 1 & 2\\ 1 & 3 \end{array}\right] \end{align} \]
[,1]
[1,] -7.333333
[2,] 10.000000 [,1]
[1,] -7.333333
[2,] 10.000000
Scilab 2023.1.0 (May 23 2023, 09:23:00)
Y =
5.
8.
25.
x =
1. 1.
1. 2.
1. 3.
ans =
-7.3333333
10.
ans =
-7.3333333
10.000000Modelo quadrático
\[ Y_i=\beta_0+\beta_1t_i+\beta_2t_{i}^{2}+\varepsilon_i\\ i=1,2,3 \]
Portanto, o modelo linear estimado por OLS é:
\[ y=16-18t+7t^2\\ t\in[1,3] \]
Note que \(t\) e \(y\) são variáveis contínuas (reais) no domínio observado.
LeastSquares[{{1, 1, 1}, {1, 2, 4}, {1, 3, 9}}, {5, 8, 25}]
quadratic fit {{1,5},{2,8},{3,25}}
Y <- c(5,
8,
25)
x <- matrix(c(1, 1, 1,
1, 2, 4,
1, 3, 9),
ncol=3,
byrow=TRUE)
solve(t(x)%*%x)%*%t(x)%*%Y [,1]
[1,] 16
[2,] -18
[3,] 7 [,1]
[1,] 16
[2,] -18
[3,] 7
Scilab 2023.1.0 (May 23 2023, 09:23:00)
Y =
5.
8.
25.
x =
1. 1. 1.
1. 2. 4.
1. 3. 9.
ans =
16.000000
-18.000000
7.0000000
ans =
16.000000
-18.000000
7.0000000Modelo cúbico
\[ Y_i=\beta_0+\beta_1t_i+\beta_2t_{i}^{2}+\beta_3t_{i}^{3}+\varepsilon_i\\ i=1,2,3 \]
Portanto, o modelo linear estimado por OLS é:
\[ y=8.97-5.11t-0.03t^2+1.17t^3\\ t\in[1,3] \]
Fit[{5, 8, 25}, {1, t, t^2, t^3}, t]
X <- c(1, 2, 3)
Y <- c(5, 8, 25)
plot(X,Y, ylim=c(min(Y)*0.5,max(Y)*1.1), xlim=c(min(X)*0.9,max(X)*1.1),
pch = 19, main = "Comparação entre Ajustes\nLinear, Quadrático e Cúbico")
c <- pracma::polyfit(X, Y, 1)
print(s <- pracma::poly2str(c))[1] "10*x^1 - 7.333333"curve(c[1]*x+c[2], min(X), max(X),
col = "black", lwd = 2, lty = 1, add=TRUE)
c <- pracma::polyfit(X, Y, 2)
print(s <- pracma::poly2str(c))[1] "7*x^2 - 18*x + 16"curve(c[1]*x^2+c[2]*x+c[3], min(X), max(X),
col = "black", lwd = 2, lty = 2, add=TRUE)
c <- pracma::polyfit(X, Y, 3)
print(s <- pracma::poly2str(c))[1] "2.666667*x^3 - 9*x^2 + 11.33333*x"curve(c[1]*x^3+c[2]*x^2 + c[3]*x+c[4], min(X), max(X),
col = "black", lwd = 2, lty = 3, add=TRUE)
# Adicionando a legenda
legend("topleft", legend = c("Dados", "Linear", "Quadrática", "Cúbica"),
col = c("black", "black", "black"), pch = c(19, NA, NA, NA),
lty = 1:3, bty = "n")
Com 3 pontos observados não é possível estimar um modelo polinomial de ordem 4.
Um critério para avaliação da curva ajustada aos pontos observados é o índice de qualidade de ajuste denominado coeficiente de determinação:
\[ R^2=\dfrac{\sum\left(y_i-\bar{Y}\right)^2}{\sum\left(Y_i-\bar{Y}\right)^2} \]
Note que \(R^2\in[0,1]\).
Quanto mais próximo de 1 o índice de qualidade de ajuste estiver, melhor o ajuste da curva aos pontos observados. Se o índice de qualidade de ajuste é exatamente igual a 1, então ocorre a interpolação.
| Modelo | \(R^2\) | \(k\) parâmetros | \(R^2/k\) |
|---|---|---|---|
| linear | 0.860 | 2 | 0.43 |
| quadrático | 1.000 | 3 | 0.33 |
| cúbico | 1.000 | 4 | 0.25 |
O critério de qualidade de ajuste relativo para comparar os modelo aos pontos observados é \(R^2/k\). Quando maior seu valor, melhor é o modelo.
Portanto, o modelo linear é o que melhor se ajusta aos pontos, pois tem o maior valor de \(R^2/k=0.43\).
O sobreajustamento (overfitting) ocorreu nos modelos quadrático e cúbico pelo fato de \(R^2=1\), i.e., o modelo “decorou” os pontos observados. Pelo princípio da parcimônia e também pela experiência prática, os modelos até a ordem 3, inclusive, costumam ser úteis nas aplicações com pelo menos quatro pontos.
Equação diferencial parcial (EDP)
Uma equação diferencial parcial (EDP) é uma equação que envolve derivadas parciais de uma função desconhecida em relação a duas ou mais variáveis independentes. Enquanto as equações diferenciais ordinárias (ODE) envolvem apenas derivadas em relação a uma única variável independente, as EDP incluem derivadas parciais em relação a várias variáveis independentes.
As EDP são amplamente utilizadas para modelar fenômenos físicos e processos dinâmicos que envolvem várias variáveis independentes. Elas são especialmente relevantes em áreas como física, engenharia, finanças, meteorologia, entre outras. As EDP descrevem relações entre a função desconhecida e suas derivadas parciais em termos das variáveis independentes e possíveis constantes ou parâmetros.
As EDP podem ser classificadas em diferentes tipos, dependendo de suas características. Alguns exemplos comuns de EDP incluem a equação da difusão, a equação da onda, a equação de Laplace e a equação de Schrödinger, entre outras. Cada tipo de EDP possui propriedades e métodos de solução específicos.
A solução de uma EDP envolve encontrar uma função que satisfaça a equação diferencial e quaisquer condições de contorno ou condições iniciais especificadas. Resolver EDP pode ser um desafio complexo e requer técnicas matemáticas avançadas, como transformadas de Fourier, métodos de separação de variáveis, métodos numéricos, entre outros.
Como as populações de moléculas, células ou organismos raramente se distribuem uniformemente em um ambiente homogêneo, seus movimentos, migrações e redistribuições são de grande interesse. No nível individual, o movimento pode resultar de processos mecanobioquímicos específicos, de contrações macroscópicas de músculos ou de fluxos ameboides. No nível populacional, esses diferentes mecanismos podem ter menos influência sobre a migração líquida do que outros fatores, como: (1) variações no ambiente, (2) densidade populacional e grau de superlotação e (3) movimento do fluido ou do ar no qual os organismos vivem.
No nível coletivo, muitas vezes é apropriado fazer uma suposição de continuidade, ou seja, representar células ou organismos discretos por meio de distribuições contínuas de densidade. Isso leva a modelos de equações diferenciais parciais (EDPs) que frequentemente são análogos aos modelos clássicos de difusão molecular, convecção ou atração.
Historicamente, os modelos biológicos envolvendo equações diferenciais parciais (EDPs) remontam aos trabalhos de K. Pearson e J. Blakeman no início do século XX. Na década de 1930, outros pesquisadores, incluindo R. A. Fisher, aplicaram EDPs para modelar a propagação espacial de genes e doenças. A década de 1950 presenciou vários desenvolvimentos importantes, incluindo o trabalho de A. M. Turing sobre formação de padrões (ver Capítulo 11) e a análise de Skellam (1951), que foi um dos primeiros a aplicar formalmente a equação de difusão na modelagem da dispersão aleatória de uma população na natureza. Alguns desses modelos, bem como outros provenientes da biologia molecular, celular e populacional, são apresentados neste capítulo.
Exemplo: Equação da onda unidimensional de d’Alembert
O primeiro artigo científico que resolveu uma equação diferencial parcial é atribuído a Jean le Rond d’Alembert. Em 1746, d’Alembert publicou um artigo chamado “Recherches sur la courbe que forme une corde tendue mise en vibration” (“Investigações sobre a curva formada por uma corda esticada posta em vibração”, em tradução livre). Nesse artigo, d’Alembert resolveu a equação diferencial parcial conhecida como a equação da onda unidimensional.
Wave equation in one space dimension: Wikipedia
A equação diferencial parcial que descreve a propagação de ondas em uma corda vibrante é dada por:
\[ \dfrac{\partial^2 u}{\partial t^2} = c^2 \dfrac{\partial^2 u}{\partial x^2} \]
Sendo que:
- \(u(x, t)\) é a função que descreve o deslocamento da onda no ponto \(x\) e no tempo \(t\).
- \(c\) é a velocidade de propagação da onda.
- \(\dfrac{\partial^2 u}{\partial t^2}\) é a aceleração da onda.
- \(\dfrac{\partial^2 u}{\partial x^2}\) representa a curvatura espacial da onda.
d’Alembert descobriu que a solução geral da equação da onda pode ser expressa como a superposição de duas funções que se propagam em direções opostas:
\[ u(x, t) = \dfrac{f(x - ct) + g(x + ct)}{2} \]
- \(f(x - ct)\) representa uma onda que se propaga para a direita com velocidade \(c\).
- \(g(x + ct)\) representa uma onda que se propaga para a esquerda com velocidade \(c\).
Essa solução mostra que a forma da onda se preserva enquanto se desloca ao longo do eixo \(x\) com velocidade \(c\).
A solução de d’Alembert reflete o princípio da superposição de ondas: a solução geral é a soma de duas ondas que se movem em direções opostas. Isso explica fenômenos como a propagação de vibrações em uma corda esticada ou ondas sonoras em um meio homogêneo.
D’Alembert obteve a solução particular da equação da onda aplicando o método de separação de variáveis. Ele mostrou que a forma geral da solução envolve a soma de duas funções, uma avançando e outra retrocedendo na direção da propagação da onda.
O trabalho de d’Alembert foi um marco importante no campo das equações diferenciais parciais e contribuiu para o desenvolvimento de métodos de solução e compreensão dos fenômenos ondulatórios. Desde então, as equações diferenciais parciais têm sido objeto de estudo e pesquisa em várias áreas da ciência e da matemática aplicada.
Para uma onda inicial dada por uma função senoidal:
\[ \begin{align} u(x, 0) &= \sin(x)\\ \dfrac{\partial u}{\partial t}(x,0)&=0 \end{align} \]
A segunda condição inicial garante que a onda começa sem velocidade inicial (onda estacionária).
A solução de d’Alembert é:
\[ u(x, t) = \dfrac{\sin(x - ct) + \sin(x + ct)}{2}=\cos(ct)\sin(x) \]
(sin(x - ct) + sin(x + ct))/2= cos(ct) sin(x)
Isso corresponde a uma onda estacionária se as duas ondas se propagam com a mesma amplitude e frequência.
Mathematica:PDE_wave.nb
DSolveValue[{D[u[x, t], {t, 2}] == c^2 D[u[x, t], {x, 2}]
&& u[x, 0] == Sin[x] && Derivative[0, 1][u][x, 0] == 0},
u[x, t], {x, t}] // Simplify
# Definindo a função da onda estacionária
onda <- function(x, t, c) {
(sin(x - c * t) + sin(x + c * t))/2
}
# Parâmetros da onda
c <- 1 # Velocidade da onda
x <- seq(0, 2 * pi, length.out = 300) # Eixo x
# Salvar como GIF
animation::saveGIF({
for (t in seq(0, 2 * pi, length.out = 50)) {
# Calcular a onda para o tempo t
y <- onda(x, t, c)
# Plotar a onda
plot(x, y, type = "l", col = "blue", lwd = 2,
ylim = c(-2, 2), xlab = "x", ylab = "u(x, t)",
main = paste("Onda Estacionária: t =", round(t, 2)))
grid()
}
}, movie.name = "onda_estacionaria.gif", interval = 0.1)Output at: onda_estacionaria.gif[1] TRUEA solução desta EDP para outras condições de contorno está em Barreira (2014, p. 520, 525-7):
DSolveValue[{D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}]
&& u[x, 0] == 0 && u[0, t] == Cos[Pi t] &&
Derivative[0, 1][u][x, 0] == 1}, u[x, t], {x, t}] // Simplify //
Simplify
Equação diferencial parcial linear
Uma equação diferencial parcial linear (EDP linear) é uma equação na qual a função desconhecida e suas derivadas parciais aparecem de forma linear, ou seja, não são elevadas a potências diferentes de 1 nem multiplicadas entre si.
Em outras palavras, se a função \(u(x, y, \dots)\) e suas derivadas parciais aparecem apenas somadas ou multiplicadas por funções conhecidas (mas não entre si), a EDP é linear.
Exemplo simples:
\[ a(x, y) \dfrac{\partial u}{\partial x} + b(x, y) \dfrac{\partial u}{\partial y} + c(x, y) \cdot u = f(x, y) \]
Aqui:
- \(a(x, y)\), \(b(x, y)\) e \(c(x, y)\) são funções conhecidas.
- \(f(x, y)\) é o termo não homogêneo.
- A função \(u\) e suas derivadas aparecem de forma linear.
A seguir há dois exemplos clássicos de equações diferenciais parciais lineares (PDEs lineares):
- Equação da Onda:
A equação da onda descreve a propagação de ondas, como ondas sonoras ou ondas em uma corda vibrante:
\[ \dfrac{\partial^2 u}{\partial t^2} = c^2 \dfrac{\partial^2 u}{\partial x^2} \]
Nesta equação:
- \(u = u(x, t)\) é a função desconhecida, representando a amplitude da onda.
- \(c\) é a velocidade de propagação da onda.
- A equação é linear porque \(u\) e suas derivadas aparecem apenas na primeira potência e não se multiplicam.
- Equação de Laplace:
A equação de Laplace aparece em problemas de equilíbrio térmico, potencial eletrostático e fluxo de fluidos:
\[ \dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2} = 0 \]
Nesta equação:
- \(u = u(x, y)\) é a função desconhecida.
- A soma das segundas derivadas parciais é igual a zero.
- A equação é linear porque não há produtos ou potências das derivadas de \(u\).
Equação de difusão unidimensional
A equação de difusão (dispersão) unidimensional, também conhecida como equação do calor, foi formulada por Joseph Fourier no início do século XIX. Fourier apresentou suas ideias sobre a propagação do calor e a equação de difusão em seu trabalho seminal “Théorie analytique de la chaleur” (Teoria Analítica do Calor), publicado em 1822. Nesse trabalho, Fourier desenvolveu a teoria matemática da condução de calor e introduziu a equação diferencial parcial que descreve a difusão de calor em sólidos.
Existem muitas combinações possíveis de substâncias e solventes que poderiam ser usadas em um experimento de injeção em um tubo cheio de líquido. Aqui estão alguns exemplos comuns:
- Azul de Metileno em Água: Usado frequentemente para estudos de difusão e transporte de massa.
- Cloreto de Sódio em Água: Comum em experimentos de química.
- Ácido Acético em Água: Usado em muitas reações químicas e estudos de difusão.
- Oxigênio em Água: Importante em estudos de transferência de gás e processos biológicos.
A escolha da substância e do solvente dependerá dos objetivos específicos do experimento e das propriedades físicas e químicas desejadas.
Suponha que injetamos uma substância em um tubo cheio de um líquido solvente (Fig. 12.5). As moléculas da substância estão em movimento aleatório. Eles se espalharão em todas as direções. Assumimos que o líquido solvente não está em movimento.
Em locais de alta concentração, as moléculas tenderão a diminuir em número e, inversamente, em locais de baixa concentração, elas tenderão a aumentar em número. Esse processo de migração aleatória, chamado difusão, finalmente terminará com as moléculas em densidades iguais em todo o tubo.
Fig. 12.5. As moléculas de uma substância dissolvida em um líquido estão em movimento aleatório. Isso resulta em uma migração lenta, chamada difusão. A densidade C(x, t) da substância dissolvida tende à uniformidade no espaço.
Na Fig. 12.5 colocamos um eixo x paralelo ao eixo do tubo. Para simplificar, assumimos que a concentração da substância varia apenas na direção x. A cada x, a concentração depende também do instante de tempo t. Assim, podemos denotar a concentração por:
\[ C=C(x,t) \]
Como unidade de medida, tomamos, por exemplo, g/cm3.
A função \(C\) é plotada na Fig. 12.5 para três diferentes instantes de tempo \(t_1\), \(t_2\), \(t_3\).
Nosso objetivo é determinar a função desconhecida \(C\) ou, pelo menos, encontrar uma equação que \(C\) deva satisfazer.
Para isso, assumimos que \(C\) é diferenciável em relação a \(x\) e a \(t\).
Estritamente falando, a matéria é discreta. No entanto, o número de moléculas é tão grande que o erro no uso de uma função suave é desprezível (cf. Fig. 9.5).
Faça \(x = 0\) coincidir com a extremidade esquerda do tubo e \(x = a\) com a extremidade direita. Denotamos a área (constante) da seção transversal por \(A\).
Seja \(M = M(x, t)\) a massa da substância contida entre a extremidade esquerda do tubo e a seção transversal em \(x\) no instante \(t\) (Fig. 12.5).
Para um primeiro resultado, consideramos \(M\) apenas como uma função de \(x\) e mantemos o tempo \(t\) constante.
Seja \(x\) um valor fixo (local de aplicação da substância) e \(\Delta x = h\) um aumento de \(x\).
A massa \(M(x + h, t)\) consiste na massa \(M(x, t)\) e na massa adicional \(\Delta M\) localizada entre \(x\) e \(x + h\) (Fig. 12.6).
Por isso:
\[ \Delta M = M(x+h, t) - M(x, t) \]
O volume ocupado por essa massa é \(A h\).
Fig. 12.6. Explicação de ΔM.
Sua densidade média é, portanto:
\[ \dfrac{\Delta M}{A h}=\dfrac{1}{A}\dfrac{M(x+h,t)-M(x,t)}{h} \]
Como \(h\) tende a zero, o lado esquerdo tende a \(C\) e o segundo fator no lado direito para \(\partial M/\partial x\).
Por isso,
\[ \begin{align} C&= \dfrac{1}{A}\dfrac{\partial M}{\partial x} \\ \dfrac{\partial M}{\partial x}&= A C \end{align} \]
Para obter um segundo resultado, consideramos \(M\) como uma função apenas de \(t\) e mantemos \(x\) constante.
A massa \(M(x, t)\) aumenta com o tempo se o fluxo líquido de moléculas em \(x\) for direcionado para a esquerda.
Esse caso ocorre se \(C\) aumenta com o aumento de \(x\), ou seja, se \(\partial C/\partial x > 0\) (ver Fig. 12.5).
Chamamos de \(\partial C/\partial x\) o gradiente da concentração.
É intuitivamente claro que \(M\) aumenta mais rapidamente à medida que esse gradiente se torna maior.
É razoável esperar que a taxa de aumento de \(\partial M/\partial t\) seja proporcional a \(\partial C/\partial x\), bem como à área \(A\) da seção transversal.
De fato, essa suposição é bem apoiada por fatos experimentais.
Portanto:
\[ \dfrac{\partial M}{\partial t}= D A \dfrac{\partial C}{\partial x} \]
sendo que \(D\) denota uma constante positiva conhecida como constante de difusão.
Um argumento semelhante é válido quando \(\partial C/\partial x < 0\).
Então o fluxo líquido é para a direita e \(M\) diminui em função do tempo.
A fórmula \(\dfrac{\partial M}{\partial t}= D A \dfrac{\partial C}{\partial x}\) contém duas funções desconhecidas, \(M\) e \(C\).
Como estamos interessados em \(C\), eliminaremos \(M\).
Para isso, derivamos \(\dfrac{\partial M}{\partial t}= D A \dfrac{\partial C}{\partial x}\) em relação a \(x\) e obtemos
\[ \dfrac{\partial^2 M}{\partial x \partial t}= D A \dfrac{\partial^2 C}{\partial x^2} \]
Em seguida, derivamos \(\dfrac{\partial M}{\partial x}= A C\) em relação a \(t\) e obtemos
\[ \dfrac{\partial^2 M}{\partial t \partial x}= A \dfrac{\partial C}{\partial t} \]
Note que \(\dfrac{\partial^2 M}{\partial x \partial t}=\dfrac{\partial^2 M}{\partial t \partial x}\).
Então:
\[ \dfrac{\partial C}{\partial t}= D \dfrac{\partial^2 C}{\partial x^2} \]
Esta é a famosa equação de difusão unidimensional.
Pode ser generalizado para todas as três dimensões espaciais \(x\), \(y\), \(z\).
DSolve[D[f[x,t],{t,1}] - D D[f[x,t],{x,2}]==0, f[x,t], {x,t}]
Exemplo 12.4.1: Primeira solução particular
Encontre os coeficientes \(a\) e \(b\) tais que
\[ C(x, t) = \exp\left({ax + bt}\right) \]
satisfaz a equação de difusão.
Aplicando a regra da cadeia obtemos:
\[ \dfrac{\partial C}{\partial t}= be^{ax + bt}\\ \dfrac{\partial^2 C}{\partial x^2}= a^2 e^{ax + bt} \]
Substituindo essas derivadas na equação de difusão obtemos:
\[ be^{ax + bt}=D a^2 e^{ax + bt} \]
Esta equação é satisfeita para todos os valores de \(x\) e \(t\) se escolhermos o coeficiente \(a\) arbitrariamente e se adotarmos \(b=Da^2\).
Portanto:
\[ C(x, t) = \exp\left({a(x + Dat)}\right) \]
Exemplo 12.4.2: Segunda solução particular
Verifique se
\[ C(x, t) = t^{-1/2} \exp\left(-\dfrac{x^2}{4Dt}\right) \]
é uma solução particular da equação de difusão.
plot (1/sqrt(t)) e^(-x^2/(4 t))
\(C(x, t)\) é um produto de duas funções de \(t\).
Assim:
\[ \dfrac{\partial C}{\partial t}=\left(t^{-1/2}\right)^{\prime} \exp\left(-\dfrac{x^2}{4Dt}\right) + t^{-1/2} \left(\exp\left(-\dfrac{x^2}{4Dt}\right)\right)^{\prime} \]
Então, obtemos:
\[ \begin{align} \dfrac{\partial C}{\partial t}&=\left(-\dfrac{1}{2}t^{-1/2}+\dfrac{x^2}{4Dt}t^{-1/2}\right) \exp\left(-\dfrac{x^2}{4Dt}\right)\\ \dfrac{\partial C}{\partial x}&=t^{-1/2}\dfrac{-x}{2Dt}\exp\left(-\dfrac{x^2}{4Dt}\right)\\ \dfrac{\partial^2 C}{\partial x^2}&=t^{-1/2}\dfrac{-1}{2Dt}\exp\left(-\dfrac{x^2}{4Dt}\right)+t^{-1/2}\left(\dfrac{-x}{2Dt}\right)^{2}\exp\left(-\dfrac{x^2}{4Dt}\right) \end{align} \]
Substituindo essas derivadas na equação de difusão obtemos:
\[ \dfrac{x^2}{4D} t^{-5/2} - \dfrac{1}{2} t^{-3/2} = D t^{-1/2} \left(\dfrac{x^2}{4 D^{2} t^{2}}-\dfrac{1}{2Dt}\right) \]
Por álgebra simples, vemos que esta equação vale para todos os valores de \(t\) e \(x\).
((x^2)/(4D)) t^(-5/2) + (-1/2) t^(-3/2) = D t^(-1/2) (((x^2)/(4 D^2 t^2))-1/(2 D t))
\(C(x, t)\) é em forma de sino para cada \(t\) fixo (cf. a fórmula para a distribuição normal na Seção 13.11). Um gráfico aproximado de \(C(x, t)\) é mostrado na Fig. 12.5.
Mesmo sem uma solução explícita, a equação de difusão revela algumas propriedades principais da solução.
Considere a Fig. 12.7a. Lá assumimos que \(C(x, t)\) é uma função linear de \(x\) entre dois pontos com abcissas \(x_1\) e \(x_2\). Isso significa que \(\partial^2 C/\partial x^2=0\). A equação de difusão então implica \(\partial C/\partial t=0\), isto é, \(C(x, t)\) não muda em função de \(t\). Portanto, temos um fluxo estacionário de moléculas para todo \(x\) em \([x_1,x_2]\).
Fig. 12.7. Três tipos diferentes de função de densidade C(x, t). Em (a) a densidade é uma função linear de x, em (b) a densidade é convexa para baixo e em (c) convexa para cima. O fluxo de moléculas depende essencialmente desta forma.
Se, no entanto, \(C(x, t)\) for côncava como na Fig. 12.7b, sabemos que a segunda derivada \(\partial^2 C/\partial x^2>0\). Em virtude da equação de difusão a taxa de aumento \(\partial C/\partial t>0\).
Portanto, \(C(x, t)\) aumentará em função de \(t\) para todo \(x\) em \([x_1,x_2]\). O resultado é bastante plausível, pois a Fig. 12.7b indica um “vale” de concentração entre \(x_1\) e \(x_2\).
Exatamente o oposto é verdadeiro na Fig. 12.7c, onde temos um “pico” de concentração. A função \(C(x, t)\) é convexa. Por isso, \(\partial^2 C/\partial x^2<0\) e \(\partial C/\partial t<0\). A densidade \(C(x, t)\) diminuirá em função do tempo para todo \(x\) em \([x_1,x_2]\).
Dentro das dimensões de uma célula viva, o transporte de massa por difusão é bastante eficiente. Atua em segundos ou no máximo em minutos. Em um corpo maior, as moléculas não podem se mover rápido o suficiente por difusão.
A distribuição de nutrientes seria muito lenta. No ser humano levaria uma vida inteira para obter o açúcar que foi alimentado no estômago para se difundir nos pés e nas mãos (Went, 1968). Portanto, é necessário um sistema de distribuição mais rápido, bem diferente da difusão. Esta é o da convecção pela corrente sanguínea. A convecção pela corrente sanguínea refere-se ao transporte de substâncias dissolvidas, como nutrientes, gases e medicamentos, pelo fluxo sanguíneo no sistema circulatório. Na prática, o estudo da convecção na corrente sanguínea é essencial para entender a distribuição de medicamentos, a troca de gases nos pulmões e o fornecimento de nutrientes às células do corpo.
\[\Diamond\]
Equação do cabo
A importância da equação diferencial de difusão vai muito além da difusão. Em ecologia, a invasão de uma grande área por uma nova espécie pode ser matematicamente tratada pelo mesmo método que o movimento aleatório de moléculas (ver Pielou, 1969, Cap. 11). Da mesma forma, nas epidemias, a disseminação de uma doença infecciosa em uma grande área segue o mesmo padrão. Uma troca de energia cinética entre moléculas vizinhas conhecida como condução de calor é tratada matematicamente da mesma forma que a difusão. Aplicações mais sofisticadas da equação de difusão são feitas em genética de populações e evolução.
Outra equação diferencial parcial biologicamente importante é a equação do cabo ou telegráfica:
\[ \dfrac{\partial^2 V}{\partial x^2}-a^2 \dfrac{\partial^2 V}{\partial t^2}+b\dfrac{\partial V}{\partial t}+cV=d \]
Aqui, \(V = V(X, t)\) denota a tensão na abcissa \(x\) do cabo e no instante de tempo \(t\).
A equação do cabo foi fundamental para a análise do primeiro cabo telegráfico transatlântico. Este cabo, lançado em 1858, permitiu a comunicação telegráfica entre a América do Norte e a Europa. Devido aos efeitos resistivos e capacitivos do cabo, o sinal elétrico se atenuava e distorcia ao longo da distância, exigindo um modelo matemático para entender e melhorar a transmissão dos sinais.
- Termo \(\dfrac{\partial^2
V}{\partial x^2}\):
- Este termo representa a difusão espacial do potencial \(V\) ao longo do cabo. É análogo ao termo de difusão em equações de condução de calor ou de difusão de partículas.
- Termo \(- a^2 \dfrac{\partial^2
V}{\partial t^2}\):
- Este termo representa a aceleração temporal do potencial \(V\). O coeficiente \(a\) está relacionado à velocidade de propagação do sinal no cabo.
- Termo \(b \dfrac{\partial
V}{\partial t}\):
- Este termo representa um amortecimento ou dissipação do sinal no tempo. O coeficiente \(b\) está associado à resistência do material do cabo.
- Termo \(c V\):
- Este termo representa um ganho ou perda proporcional ao próprio potencial \(V\). Pode estar relacionado a fontes ou perdas no sistema.
- Termo \(d\):
- Este termo é uma fonte ou força externa que não depende nem de \(x\) nem de \(t\).
Hoje, a equação do cabo é aplicada na teoria da condução nervosa (ver Beier, 1962, Cole, 1968, Yihe & Timofeeva, 2020).
Sempre que flutuações aleatórias estão sendo consideradas no tratamento de processos biológicos, a teoria da probabilidade entra em cena.
O método matemático apropriado é o uso de processos estocásticos.
Eles ocorrem em trabalhos teóricos sobre crescimento populacional, competição entre espécies, disseminação de infecções etc.
Entre as ferramentas matemáticas necessárias para o estudo de processos estocásticos está uma grande variedade de equações diferenciais parciais. Aqui temos que pular detalhes. O leitor é remetido a Bailey (1964) e a Chiang (1968).
Finanças Quantitativas: Modelo de Black-Scholes-Merton
A equação de Black-Scholes-Merton é uma equação diferencial parcial (EDP) fundamental na matemática financeira para a modelagem do preço de opções.
Em 1997, o Prêmio Nobel de Ciências Econômicas foi concedido a Myron Scholes e Robert Merton pelo seu trabalho no desenvolvimento do modelo de precificação de opções, que inclui a equação de Black-Scholes-Merton. Fischer Black, que também contribuiu significativamente para o desenvolvimento da equação, Ele faleceu em 1995 e, portanto, não foi elegível para o prêmio. O trabalho deles revolucionou os mercados financeiros, permitindo uma precificação mais precisa dos derivativos e promovendo uma maior compreensão do risco financeiro.
- Opção de Compra (Call Option)
A opção de compra \(C\) pode ser definida como:
\[ C_T = \max(S_T - K, 0) \]
Onde:
\(C_T\) = valor da opção de compra na data de vencimento
\(S_T\) = preço do ativo subjacente na data de vencimento
\(K\) = preço de exercício (strike price)
Opção de Venda (Put Option)
A opção de venda \(P\) pode ser definida como:
\[ P_T = \max(K - S_T, 0) \]
Onde:
- \(P_T\) = valor da opção de venda na data de vencimento
- \(S_T\) = preço do ativo subjacente na data de vencimento
- \(K\) = preço de exercício (strike price)
Estas fórmulas representam o valor das opções no vencimento, onde o valor é o maior entre a diferença entre o preço do ativo subjacente e o preço de exercício, ou zero.
A equação de Black-Scholes-Merton para o preço racional de uma opção \(V(S_t, t)\) é dada por:
\[ \dfrac{\partial V}{\partial t} + \dfrac{1}{2} \sigma^2 S^2 \dfrac{\partial^2 V}{\partial S^2} + r S \dfrac{\partial V}{\partial S} - r V = 0 \]
onde:
- \(V(S, t)\) é o preço da opção como função do preço do ativo subjacente \(S\) e do tempo \(t\).
- \(\sigma\) é a volatilidade do ativo subjacente.
- \(r\) é a taxa de juros livre de risco.
- \(\dfrac{\partial V}{\partial t}\) é a derivada parcial de \(V\) em relação ao tempo \(t\).
- \(\dfrac{\partial^2 V}{\partial S^2}\) é a derivada parcial de segunda ordem de \(V\) em relação ao preço do ativo subjacente \(S\).
- \(\dfrac{\partial V}{\partial S}\) é a derivada parcial de primeira ordem de \(V\) em relação ao preço do ativo subjacente \(S\).
Claro, aqui estão as fórmulas da solução para a equação de Black-Scholes-Merton:
Para uma opção de compra (call):
\[ C(S_t, t) = S_t N(d_1) - K e^{-r(T-t)} N(d_2) \]
Para uma opção de venda (put):
\[ P(S_t, t) = K e^{-r(T-t)} N(-d_2) - S_t N(-d_1) \]
Onde:
\[ \begin{align} d_1 &= \dfrac{\ln(S_t / K) + \left(r + \dfrac{\sigma^2}{2}\right)(T - t)}{\sigma \sqrt{T - t}}\\ d_2 &= d_1 - \sigma \sqrt{T - t} \end{align} \]
- \(S_t\) é o preço atual do ativo subjacente.
- \(K\) é o preço de exercício da opção.
- \(r\) é a taxa de juros livre de risco.
- \(\sigma\) é a volatilidade do ativo subjacente.
- \(T - t\) é o tempo até a maturidade da opção.
- \(N(\cdot)\) é a função de distribuição cumulativa da distribuição normal padrão.
PDEs importantes
As equações diferenciais parciais (PDEs) são fundamentais na descrição de diversos fenômenos naturais e tecnológicos. Cada PDE clássica tem sua origem e evolução histórica, com soluções que foram desenvolvidas ao longo dos séculos.
1. Equação de Laplace
Introdução Histórica: A equação de Laplace é nomeada em homenagem ao matemático francês Pierre-Simon Laplace, que a introduziu no final do século XVIII. Ela surgiu no contexto da teoria do potencial e é fundamental em eletrostática, gravitação e teoria do potencial.
Aplicações:
- Física: Potenciais gravitacionais e eletrostáticos.
- Química: Difusão de substâncias em equilíbrio.
- Biologia: Distribuição de temperatura em organismos e difusão de moléculas em células.
Equação:
\[ \nabla^2 u = 0 \quad \text{ou} \quad \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0 \]
Solução:
Para uma função \(u(x, y)\) em duas dimensões:
\[ u(x, y) = \sum_{n=1}^{\infty} \left( A_n e^{n \pi x} + B_n e^{-n \pi x} \right) \left( C_n \sin(n \pi y) + D_n \cos(n \pi y) \right) \]
2. Equação de Poisson
Introdução Histórica: A equação de Poisson é uma extensão da equação de Laplace, introduzida por Siméon Denis Poisson no início do século XIX. Ela é essencial na eletrostática e em problemas de potencial.
Aplicações:
- Física: Eletrostática e teoria do potencial.
- Química: Distribuição de carga em moléculas.
- Biologia: Modelagem de campos de concentração de nutrientes em tecidos.
Equação:
\[ \nabla^2 u = f(x, y) \quad \text{ou} \quad \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x, y) \]
Solução:
Para \(f(x, y) = k\), onde \(k\) é constante:
\[ u(x, y) = \frac{k}{2}(x^2 + y^2) + Ax + By + C \]
3. Equação da Onda
Introdução Histórica: A equação da onda foi formulada no século XVIII por Jean le Rond d’Alembert, que a utilizou para descrever a propagação de ondas em cordas vibrantes.
Aplicações:
- Física: Propagação de ondas sonoras, sísmicas e eletromagnéticas.
- Química: Dinâmica de reações em meios oscilantes.
- Biologia: Propagação de sinais em nervos e ondas de pressão em fluidos biológicos.
Equação:
\[ \frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u \quad \text{ou} \quad \frac{\partial^2 u}{\partial t^2} = c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \]
Solução:
Para uma onda unidimensional em uma corda com condições de contorno \(u(0, t) = u(L, t) = 0\):
\[ u(x, t) = \sum_{n=1}^{\infty} \left( A_n \cos \left( \frac{n \pi c t}{L} \right) + B_n \sin \left( \frac{n \pi c t}{L} \right) \right) \sin \left( \frac{n \pi x}{L} \right) \]
4. Equação do Calor
Introdução Histórica: A equação do calor foi introduzida por Joseph Fourier no início do século XIX, durante seus estudos sobre a condução de calor.
Aplicações:
- Física: Condução de calor em sólidos.
- Química: Difusão de substâncias em meios homogêneos.
- Biologia: Distribuição de temperatura em organismos vivos e difusão de nutrientes em células.
Equação:
\[ \frac{\partial u}{\partial t} = \alpha \nabla^2 u \quad \text{ou} \quad \frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \]
Solução: Para uma barra unidimensional de comprimento \(L\) com condições de contorno \(u(0, t) = u(L, t) = 0\):
\[ u(x, t) = \sum_{n=1}^{\infty} B_n \sin \left( \frac{n \pi x}{L} \right) e^{-\left( \frac{n \pi}{L} \right)^2 \alpha t} \]
5. Equação de Schrödinger
Introdução Histórica: A equação de Schrödinger foi formulada por Erwin Schrödinger em 1925 e é um pilar da mecânica quântica, descrevendo a evolução temporal de estados quânticos.
Aplicações:
- Física: Comportamento de partículas quânticas.
- Química: Estrutura eletrônica de átomos e moléculas.
- Biologia: Fenômenos quânticos em processos biológicos como fotossíntese e visão.
Equação:
\[ i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + V\psi \]
Solução:
Para uma partícula em um poço de potencial infinito unidimensional:
\[ \psi_n(x, t) = \sqrt{\frac{2}{L}} \sin \left( \frac{n \pi x}{L} \right) e^{-i E_n t / \hbar} \] onde \(E_n = \frac{n^2 \pi^2 \hbar^2}{2mL^2}\).
6. Equação de Klein-Gordon
Introdução Histórica: A equação de Klein-Gordon, formulada por Oskar Klein e Walter Gordon em 1926, é uma extensão relativística da equação de Schrödinger.
Aplicações:
- Física: Descrição de partículas relativísticas.
- Química: Fenômenos quânticos em altas energias.
- Biologia: Fenômenos relativísticos em biologia teórica.
Equação:
\[ \frac{\partial^2 \phi}{\partial t^2} - c^2 \nabla^2 \phi + \frac{m^2 c^4}{\hbar^2} \phi = 0 \]
Solução:
Para uma solução em uma dimensão espacial:
\[ \phi(x, t) = e^{i(kx - \omega t)} \]
onde \(\omega = \sqrt{c^2 k^2 + \frac{m^2 c^4}{\hbar^2}}\).
Essas são algumas das PDEs mais importantes, com suas formas explícitas, soluções analíticas e aplicações em biologia, química e física. Cada uma delas desempenha um papel crucial na compreensão e modelagem de fenômenos complexos em diversas áreas da ciência.
PDES importantes: Soluções em Mathematica
1. Equação de Laplace
Equação:
\[ \nabla^2 u = 0 \quad \text{ou} \quad \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0 \]
Comando em Mathematica:
DSolve[{Laplacian[u[x, y], {x, y}] == 0}, u[x, y], {x, y}]
Resultado:
\[ u(x, y) = C[1](I x + y) + C[2](-I x + y) \]
2. Equação de Poisson
Equação:
\[ \nabla^2 u = k \quad \text{onde} \quad k \text{ é constante} \]
Comando em Mathematica:
DSolve[{Laplacian[u[x, y], {x, y}] == k}, u[x, y], {x, y}]
Resultado:
\[ u(x, y) = -\frac{k}{4} (x^2 + y^2) + C[1] x + C[2] y + C[3] \]
3. Equação da Onda
Equação:
\[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]
Comando em Mathematica:
DSolve[{D[u[x, t], {t, 2}] == c^2 D[u[x, t], {x, 2}]}, u[x, t], {x, t}]
Resultado:
\[ u(x, t) = C[1]\left(t - \frac{x}{\sqrt{c^2}}\right) + C[2]\left(t + \frac{x}{\sqrt{c^2}}\right) \]
4. Equação do Calor
Equação: \[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]
Comando em Mathematica:
DSolve[{D[u[x, t], t] == α D[u[x, t], {x, 2}]}, u[x, t], {x, t}]
Resultado:
\[ u(x, t) = 1 + \cosh(C[1] + x C[2] + t \alpha C[2]^2) + \sinh(C[1] + x C[2] + t \alpha C[2]^2) \]
5. Equação de Schrödinger
Equação:
\[ i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + V\psi \]
Vamos adotar \(V(x) = V_0\) constante.
Comando em Mathematica:
DSolve[{I hbar D[ψ[x, t], t] == -((hbar^2)/(2 m)) D[ψ[x, t], {x, 2}] + V0 ψ[x, t]}, ψ[x, t], {x, t}]
Resultado:
\[ ψ(x, t) = e^{-\frac{i t (V0 + \frac{\hbar^2 k^2}{2 m})}{\hbar}} (C[1] e^{i k x} + C[2] e^{-i k x}) \]
6. Equação de Klein-Gordon
Equação:
\[ \frac{\partial^2 \phi}{\partial t^2} - c^2 \frac{\partial^2 \phi}{\partial x^2} + \frac{m^2 c^4}{\hbar^2} \phi = 0 \]
Comando em Mathematica:
DSolve[{D[ϕ[x, t], {t, 2}] - c^2 D[ϕ[x, t], {x, 2}] + (m^2 c^4 / ħ^2) ϕ[x, t] == 0}, ϕ[x, t], {x, t}]
Resultado:
\[ ϕ(x, t) = e^{i(kx - \omega t)}, \quad \omega = \sqrt{c^2 k^2 + \frac{m^2 c^4}{\hbar^2}} \]
Resumo das Soluções Ajustadas
- Equação de Laplace: \(u(x, y) = C[1](I x + y) + C[2](-I x + y)\)
- Equação de Poisson: \(-\frac{k}{4} (x^2 + y^2) + C[1] x + C[2] y + C[3]\)
- Equação da Onda: \(C[1](t - \frac{x}{\sqrt{c^2}}) + C[2](t + \frac{x}{\sqrt{c^2}})\)
- Equação do Calor: \(1 + \cosh(C[1] + x C[2] + t \alpha C[2]^2) + \sinh(C[1] + x C[2] + t \alpha C[2]^2)\)
- Equação de Schrödinger: \(e^{-\frac{i t (V0 + \frac{\hbar^2 k^2}{2 m})}{\hbar}} (C[1] e^{i k x} + C[2] e^{-i k x})\)
- Equação de Klein-Gordon: \(e^{i(kx - \omega t)}, \quad \omega = \sqrt{c^2 k^2 + \frac{m^2 c^4}{\hbar^2}}\)
Matemática no YouTube
- blackpenredpen: delta y vs. dy (differential)
- Physics for Students- Unleash your power!!: What is the meaning of Greek symbols | Delta, del, d | Greek letters in mathematics | Greek symbols
- Dr. Trefor Bazett : What are the big ideas of Multivariable Calculus?? Full Course Intro
- Dr. Trefor Bazett : A full course playlist for Multivariable Calculus (aka Calculus III)
- Richard Behiel: The Beauty of Linear Regression (How to Fit a Line to your Data)
- 3Blue1Brown: Mas o que é uma equação diferencial parcial? | Visão geral de Equações Diferenciais, Capítulo 2
- 3Blue1Brown: Solving the heat equation | DE3
- 3Blue1Brown: Solving the heat equation | DE3
Referências
- Allen, LJS (2007) An Introduction to Mathematical Biology. NJ: Pearson.
- Barreira, S (2014) Matemática aplicada às ciências farmacêuticas com Excel. Vol. 2. Lisboa: Escolar Editora.
- Batschelet, E (1979) Introduction to mathematics for life scientists. 3rd ed. NY: Springer.
- Batschelet, E (1978) Introdução à matemática para biocientistas. Tradução da 2ª ed. São Paulo: EDUSP e Rio de Janeiro: Interciência.
- Beuter, A et al. (2003) Nonlinear dynamics in physiology and medicine. USA: Springer.
- Edelstein-Keshet, L (2005) Mathematical Models in Biology. USA: Society for Industrial and Applied Mathematics (SIAM).
- Garfinkel, A et al. (2017) Modeling Life: The Mathematics of Biological Systems. USA: Springer.
- Giordano, FR et al. (2015) A first course in mathematical modeling. 5th ed. OH: Thomson.
- Goldman, R (2005) Curvature formulas for implicit curves and surfaces. Computer Aided Geometric Design 22(7): 632-58. https://doi.org/10.1016/j.cagd.2005.06.005.
- González-Vega, L et al. (2023) Eigenvalues of Real Matrices with Prescribed Principal Minors Sign and Descartes Law of Signs. arXiv. https://doi.org/10.48550/arXiv.2305.08861.
- Hariki, S & Abdonour, OJ (2003) Matemática aplicada em administração, economia e contabilidade. SP: Saraiva.
- Hoffmann, L, Bradley, G, Sobecki, D, & Price, M (2013) Applied Calculus for Business, Economics, and the Social and Life Sciences: Expanded Edition (11ª ed.). USA: McGraw-Hill.
- Jones, DS, Plank, MJ & Sleeman, BD (2009) Differential Equations and Mathematical Biology. 2nd ed. UK: CRC.
- Keener, J & Sneyd, J (2009) Mathematical Physiology: Systems Physiology. 2nd ed. USA: Springer.
- Moylan, PJ (1977) Matrices with positive principal minors. Linear Algebra and its Applications 17(1): 53-8. https://doi.org/10.1016/0024-3795(77)90040-4.
- Neuhauser, C & Roper, ML (2018) Calculus for Biology and Medicine. 4th ed. USA: Pearson.
- Poularikas, AD (1999) The Handbook of Formulas and Tables for Signal Processing. Boca Raton: CRC/Springer/IEEE.
- Siqueira, JO (2012) Fundamentos de métodos quantitativos: aplicados em Administração, Economia, Contabilidade e Atuária usando WolframAlpha e SCILAB. SP: Saraiva. Soluções dos exercícios em https://www.researchgate.net/publication/326533772_Fundamentos_de_metodos_quantitativos_-_Solucoes.
- Venit, S & Katz, R (2000) Eigenvalues of Matrices of Low Rank. The College Mathematics Journal 31(3): 208–10. https://doi.org/10.2307/2687491
- Teichroew, D (1964) An introduction to management science: deterministic models. NJ: Wiley.
- Yihe, L & Timofeeva, Y (2020) Exact solutions to cable equations in branching neurons with tapering dendrites. J. Math. Neurosc. 10(1). https://doi.org/10.1186/s13408-020-0078-z
- Zwillinger, D (1997) Handbook of differential equations. 3rd ed. USA: Academic.