1.8. Simulações numéricas de modelos de crescimento natural

Hora de implementarmos os métodos de Euler descritos na seção anterior.

Vamos ver variações do modelo de crescimento natural. Primeiramente, vamos considerar o modelo clássico, de uma equação diferencial ordinária $\mathrm{d}x/\mathrm{d}t = \mu x,$ com um parâmetro $\mu$ constante, representando a taxa de crescimento específico. Depois, vamos permitir que $\mu$ seja uma variável aleatória, representando uma incerteza no parâmetro. Em seguida, vamos considerar o caso em que $\mu = \{\mu_t\}_{t\geq 0}$ é um processo aleatório, representando uma variabilidade temporal aleatória no parâmetro. Por último, vamos considerar $\mu$ novamente constante mas com a equação perturbada por um termo estocástico, $\mathrm{d}x = \mu x \mathrm{d}t + \sigma \mathrm{d}W_t.$

Usaremos apenas dois pacotes do Julia. Um é o Plots.jl, para visualização das soluções. O outro é o Random, da biblioteca padrão, apenas para fixar a reprodução das notas (gerar sempre os mesmos conjuntos de números aleatórios, por motivos didáticos e de controle de versão).

using Plots
theme(:ggplot2)
using Random
rng = Xoshiro(123)

Modelo de equação diferencial ordinária

Consideramos o problema de crescimento populacional

$$\begin{cases} \displaystyle \frac{\mathrm{d}x}{\mathrm{d}t} = \mu x, \quad t > 0,\\ x(0) = x_0. \end{cases}$$

Para cada $\mu > 0$ fixo, temos a solução $x(t) = x_0 e^{\mu t},$ para $t \geq 0.$

Resolvendo o problema de valor inicial pelo método de Euler, em um intervalo $t_0 \leq t \leq T,$ com $t_0 = 0,$ $T > 0$ e um passo de tempo $\Delta t > 0,$ obtemos a sequência de valores aproximados

$$x_{n+1} = x_n + \mu \Delta t x_n, \quad n = 0, 1, \dots, N,$$

onde $N$ é um inteiro mais próximo de $T/\Delta t.$

Escolhemos valores para $T,$ $\mu$ e $x_0,$ além de $t_0=0$:

t0 = 0.0
T = 20.0
μ = 0.1
x₀ = 0.2
nothing

Escolhemos, também, um valor para o passo de tempo, obtendo a malha temporal:

Δt = 0.2
t = t0:Δt:T

Isso nos dá $N+1$ pontos na malha, onde

N = Int(round(T/Δt))

Podemos pré-alocar um vetor para a solução numérica:

x = Vector{Float64}(undef, N+1)
nothing

Iniciando o vetor com o valor de $x_0$ e iterando de acordo com o método de Euler, obtemos a solução numérica aproximada:

x[1] = x₀
for n in 2:N+1
    x[n] = (1.0 + μ * Δt)x[n-1]
end

De posse da solução podemos visualizar o crescimento exponencial e compará-lo com a solução exata

plot(t, x, title="Crescimento exponencial \$(x_0 = $x₀, \\mu = $μ, T = $T, \\Delta t = $Δt)\$", titlefont = 10, xlabel = "t", ylabel="x", label="aproximação")

plot!(t, t -> x₀ * exp(μ * t), label="sol. exata")

Refinando a malha, obtemos uma solução mais próxima, onde quase não vemos diferença:

Δt = 0.02
t = t0:Δt:T
N = Int(round(T/Δt))
x = Vector{Float64}(undef, N+1)

x[1] = x₀
for n in 2:N+1
    x[n] = (1.0 + μ * Δt)x[n-1]
end

plot(t, x, title="Crescimento exponencial \$(x_0 = $x₀, \\mu = $μ, T = $T, \\Delta t = $Δt)\$", titlefont = 10, xlabel = "t", ylabel="x", label="aproximação", color=1)

Modelo com parâmetros aleatórios

Agora, vamos assumir uma incerteza em relação ao parâmetro, dada por uma distribuição normal, com média $\bar\mu$ e desvio padrão $\sigma$: $\mu = \mathcal{N}(\bar\mu, \sigma^2).$

Façamos uma amostragem de um certo número $M$ de valores.

M = 100
μ̄ = 0.1
σ = 0.02
μ = μ̄ .+ σ * randn(rng, M)

histogram(μ, bins = 20, xlims=(0.0, 0.2), title="Histogram das realizações de \$\\mu \\sim \\mathcal{N}(\\bar\\mu, \\sigma^2); \\bar\\mu = $μ̄, \\sigma=$σ \$", xlabel="\$\\mu\$", titlefont=10, label=false)

Agora, vamos resolver a equação para cada valor sorteado.

Δt = 0.2
t = t0:Δt:T
N = Int(round(T/Δt))
x = Array{Float64}(undef, N+1, M)

x[1, :] .= x₀
for n in 2:N+1
    x[n, :] .= (1.0 .+ μ * Δt) .* x[n-1, :]
end

plot(t, x, alpha = 0.2, title="Conjunto de soluções \$(x_0 = $x₀, \\mu \\sim \\mathcal{N}($μ̄, {$σ}^2), T = $T, \\Delta t = $Δt)\$", titlefont = 10, xlabel = "t", ylabel="x", label=permutedims(["soluções"; fill(nothing, M-1)]), color=1)
plot!(t, x[:, 1], label="uma realização", linewidth=2, color=2)
plot!(t, sum(x, dims=2)/size(x, 2), linewidth=3, color=:orange, label="média")

Modelo de equação diferencial aleatória

Consideramos, agora, a equação aleatória

$$\frac{\mathrm{d}X_t}{\mathrm{d}t} = \mu_t X_t,$$

onde $\{\mu_t\}_t$ é um processo aleatório dado por $\mu_t = \bar\mu + \sigma \sin(W_t),$ onde $\bar\mu, \sigma > 0$ e $\{W_t\}_{t\geq 0}$ é tal que $W_0 = 0,$ e $W_{t + \tau} - W_t \sim \mathcal{N}(0, \tau)$ independentes.

μ̄ = 0.1
σ = 0.05
wt = 0.0

for m in 1:M
    x[1, m] = x₀
    for n in 2:N+1
        global wt += randn(rng) * √Δt
        μt = μ̄ + σ * sin(wt)
        x[n, m] = (1.0 + μt * Δt) .* x[n-1, m]
    end
end

plot(t, x, alpha = 0.2, title="Soluções RODE \$(x_0 = $x₀, \\mu_t = $μ̄ + $σ\\sin(W_t), T = $T, \\Delta t = $Δt)\$", titlefont = 10, xlabel = "t", ylabel="x", label=permutedims(["soluções"; fill(nothing, M-1)]), color=1)
plot!(t, x[:, 1], label="uma realização", linewidth=2, color=2)
plot!(t, sum(x, dims=2)/size(x, 2), linewidth=3, color=:orange, label="média")

Modelo de equação diferencial estocástica

Finalmente, vamos considerar a equação estocástica

$$\mathrm{d}X_t = \bar\mu X_t \mathrm{d}t + \sigma X_t \mathrm{d}W_t,$$

onde $\{W_t\}_t$ é um processo aleatório. Mais especificamente, vamos assumir que $\{W_t\}_t$ modela um movimento Browniano, tendo incrementos independentes e estacionários, dados por $\Delta W_t = W_{t + \Delta t} - W_t \sim \mathcal{N}(0, \Delta t).$

μ̄ = 0.1
σ = 0.05

for m in 1:M
    x[1, m] = x₀
    for n in 2:N+1
        dw = randn(rng) * √Δt
        x[n, m] = (1.0 + μ̄ * Δt) * x[n-1, m] + σ * x[n-1, m] * dw
    end
end

plot(t, x, alpha = 0.2, title="Soluções SDE \$(x_0 = $x₀, \\bar\\mu = $μ̄, \\sigma = $σ, dW_t \\sim \\mathcal{N}(0, \\Delta t), T = $T, Δt = $Δt)\$", titlefont = 9, xlabel = "t", ylabel="x", label=permutedims(["soluções"; fill(nothing, M-1)]), color=1)
plot!(t, x[:, 1], label="uma realização", linewidth=2, color=2)
plot!(t, sum(x, dims=2)/size(x, 2), linewidth=3, color=:orange, label="média")

Exercícios

1. Resolva, via método de Euler, a equação diferencial

$$\frac{\mathrm{d}x}{\mathrm{d}t} = (\alpha - \beta x) x, \quad x(0) = x_0,$$

com $\alpha, \beta > 0, 0 < x_0 < \alpha / \beta, t_0 = 0, T > 0, 0 < \Delta t \ll T$ de sua escolha. Trace o gráfico da solução. Escolha parâmetros apropriados para que a mudança de concavidade, da solução, seja visível.

2. Resolva, via método de Euler, a equação diferencial

$$\mathrm{d}X_t/\mathrm{d}t = (\alpha - \beta X_t) X_t, \quad X_0 = \mathcal{N}(\bar{x}_0, \sigma_{x_0}^2),$$

onde

$$\alpha \sim \mathcal{N}(\bar\alpha, \sigma_\alpha^2), \beta \sim \mathcal{N}(\bar\beta, \sigma_\beta^2),$$

com $\bar\alpha, \sigma_\alpha, \bar\beta, \sigma_\beta, \bar{x}_0, \sigma_{x_0} > 0, t_0 = 0, T > 0, 0 < \Delta t \ll T$ de sua escolha. Trace o gráfico de um conjunto de realizações dos parâmetros. Escolha parâmetros apropriados para que boa parte das realizações exiba a mudança de concavidade.

3. Resolva, via método de Euler, a equação diferencial aleatória

$$\mathrm{d}X_t/\mathrm{d}t = (A_t - b X_t) X_t, \quad X_0 = x_0,$$

onde $A_t = a + \sigma\sin(W_t)$ e $W_{t + \tau} - W_t \sim \mathcal{N}(0, \tau),$ com $a, b, \sigma > 0,$ $0 < x_0 < a / b,$ $t_0 = 0,$ $T > 0$ e $0 < \Delta t \ll T$ de sua escolha. Trace o gráfico de um conjunto de realizações dos parâmetros. Escolha parâmetros apropriados para que a média das realizações exiba a mudança de concavidade.

4. Resolva, via método de Euler, a equação diferencial estocástica

$$\mathrm{d}X_t = (\alpha - \beta X_t) X_t \mathrm{d}t + \sigma X_t \mathrm{d}W_t, \quad X_0 = x_0,$$

onde $W_{t + \tau} - W_t \sim \mathcal{N}(0, \tau),$ com $\alpha, \beta > 0,$ $0 < x_0 < \alpha / \beta,$ $\sigma > 0,$ $t_0 = 0,$ $T > 0$ e $0 < \Delta t \ll T$ de sua escolha. Trace o gráfico de um conjunto de realizações dos parâmetros. Novamente, escolha parâmetros apropriados para que a média das realizações exiba a mudança de concavidade.