Equações de Lotka-Volterra Estocásticas

Grupo: André Luis Della Valentina, Lucas dos Santos Assmann, Vinícius Bayne Müller

Introdução

O modelo de Lotka-Volterra foi desenvolvido originalmente na década de 1920, de maneira independente por Vito Volterra e Alfred Lotka, e é utilizado para descrever a dinâmica de populações com relações de predatismo. Em sua forma mais simples, as equações de Lotka-Volterra podem ser escritas como

${\begin{cases}{\dot {x}}=x(a-by)\\{\dot {y}}=y(-c+fx)\end{cases}}$

onde $x$ e $y$ denotam, respectivamente, a densidade populacional de presas e de predadores, e $a$, $b$, $c$ e $f$ são constantes positivas.

Pode-se interpretar os parâmetros da seguinte maneira:

• $a$: taxa de crescimento livre da presa;
• $b$: taxa de predação;
• $c$: taxa de mortalidade livre do predador;
• $f$: taxa de crescimento do predador devido à predação.

É interessante notar que o sistema apresenta um ponto fixo não trivial em $(x^{\ast },y^{\ast })=\left({\frac {c}{f}},{\frac {a}{b}}\right)$. Pode-se mostrar também que as demais soluções (além da trivial) são órbitas fechadas no espaço de fase.

Nesse modelo simples, não há competição entre indivíduos de uma mesma espécie e não há limite ecológico para o sustento das populações; ou seja, a população de presas cresce exponencialmente na ausência de predadores. A fim de tornar esse modelo mais realista, pretendemos estudar versões estocásticas do mesmo, que podem ser construídas de diferentes maneiras, e analisar o efeito do ruído sobre o comportamento dinâmico.

Modelo estocástico

Buscamos construir um modelo estocástico de Lotka-Volterra utilizando um processo de Wiener; contudo, a adição de ruído branco nas equações diferenciais pode ser feita de várias maneiras. Consideraremos, aqui, duas delas: ruído adicionado nos parâmetros e ruído externo ao sistema.

Ruído externo

Como primeiro modelo, simplesmente adicionaremos um ruído externo nas equações diferenciais. Esse ruído pode ser interpretado como fatores ambientais independentes das populações, causando tanto benefícios (abundância de alimento, condições climáticas favoráveis para reprodução etc.) como prejuízos (escassez, condições climáticas desfavoráveis etc.). Como as equações tratam de densidades populacionais, devemos utilizar um ruído _multiplicativo_, pois os efeitos dos fatores externos sentidos pelas populações devem ser proporcionais ao tamanho dela. Assim, escrevemos o modelo como

${\begin{cases}{\dot {X}}=X(a-bY)+\epsilon _{x}\xi _{x}X\\{\dot {Y}}=Y(-c+fX)+\epsilon _{y}\xi _{y}Y\end{cases}}$

onde os $\xi _{i}$ são os ruídos brancos de cada variável, e $\epsilon _{i}$, a intensidade dos mesmos.

Itô

Executando um processo análogo ao realizado na seção seguinte obtém-se as equações diferenciais estocásticas no sentido de Itô:

${\begin{cases}dX=X(a-bY+{\frac {1}{2}}{\epsilon _{x}}^{2})dt+\epsilon _{x}X\bullet dW_{x}\\dY=Y(-c+fX+{\frac {1}{2}}{\epsilon _{y}}^{2})dt+\epsilon _{y}Y\bullet dW_{y}\\\end{cases}}$

Observa-se que, além de um incremento de Wiener proporcional a $X$ e $Y$, obtém-se parâmetros determinísticos "efetivos" $a\to a+{\frac {1}{2}}{\epsilon _{x}}^{2}$ e $c\to c-{\frac {1}{2}}{\epsilon _{y}}^{2}$.

Ruído nos parâmetros

Uma outra maneira de explorar o ruído nas equações de Lotka-Volterra é substituir os parâmetros determinísticos do modelo por parâmetros estocásticos. Analisaremos esse modelo proposto no sentido de Itô e no sentido de Stratonovich.

Itô

Para este caso, colocaremos o ruído de tal maneira que $a\to a+{\epsilon _{x}}X\xi _{x}(t)$ e $c\to c+{\epsilon _{y}}Y\xi _{y}(t)$, onde $\epsilon _{x}$ e $\epsilon _{y}$ são constantes e $\xi _{i}(t)$ é o ruído branco definido pelo processo de Wiener. Esse novo modelo descreveria um sistema com populações em níveis vizinhos da cadeia alimentar, abrangendo uma biodiversidade representada por parâmetros efetivos estocásticos nas taxas $a$ e $c$, proporcional ao número de presas e predadores. Ou seja, há uma diversidade de taxas de reprodutibilidade das presas e mortalidades dos predadores, o que representaria diferentes populações de diferentes espécies.

A nova equação diferencial é dada, então, por:

${\begin{cases}{\dot {X}}=X(a-bY)+\epsilon _{x}\xi _{x}(t)X^{2}\\{\dot {Y}}=Y(-c+fX)-\epsilon _{y}\xi _{y}(t)Y^{2}\end{cases}}$

Seja $A_{x}(X(t),Y(t))=X(a-bY)$, $A_{y}(X(t),Y(t))=Y(-c+fX)$, $B_{x}(X(t))=\epsilon _{x}X^{2}$ e $B_{y}(Y(t))=\epsilon _{x}Y^{2}$. Integra-se essa equação em ambos os lados de tal forma que:

${\begin{cases}\int _{t}^{t+\Delta t}{dX}=\int _{t}^{t+\Delta t}{X(a-bY)dt'}+\epsilon _{x}\int _{t}^{t+\Delta t}{X^{2}\xi _{x}(t')dt'}\\\int _{t}^{t+\Delta t}{dY}=\int _{t}^{t+\Delta t}{Y(-c+fX)dt'}+\epsilon _{y}\int _{t}^{t+\Delta t}{Y^{2}\xi _{y}(t')dt'}\end{cases}}$

A parte determinista é integrada trivialmente, porém para integrar o segundo termo das equações é necessário aproximar $B_{x}(X(t))\to B_{x}\left(X+{\frac {1}{2}}\Delta X\right)$ no intervalo $\Delta t$, e o mesmo para $B_{y}(Y(t))$. Expande-se $B_{i}$ em série de Taylor, conforme o cálculo de Itô, e então obtém-se as equações diferenciais estocásticas ao tomar o limite $\Delta t\to dt$:

${\begin{cases}dX=X(a-bY+{\epsilon _{x}}^{2}X^{2})dt+\epsilon _{x}X^{2}\bullet dW_{x}\\dY=Y(-c+fX+{\epsilon _{y}}^{2}Y^{2})dt-\epsilon _{y}Y^{2}\bullet dW_{y}\\\end{cases}}$

Ao integrá-las de $t$ a $t+\Delta t$ obtém-se o método numérico utilizado para fazer as simulações computacionais. Nota-se que a integral dos segundos termos são aproximadas para $B_{i}\Delta W_{i}$ no intervalo $\Delta t$. Então, é utilizado o método de Heun para integrar a parte determinística da equação, obtendo portanto o seguinte método:

${\begin{cases}X_{n+1}=X_{n}+{\frac {1}{2}}(k_{1x}+k_{2x})\Delta t+\epsilon _{x}X_{n}^{2}\Delta W_{x}\\Y_{n+1}=Y_{n}+{\frac {1}{2}}(k_{1y}+k_{2y})\Delta t-\epsilon _{y}Y_{n}^{2}\Delta W_{y}\\\end{cases}}$

onde

${\begin{cases}k_{1x}=X_{n}(a-bY_{n}+{\epsilon _{x}}^{2}X_{n}^{2})\\k_{2x}=(X_{n}+k_{1x}\Delta t)(a-b(Y_{n}+k_{1y}\Delta t)+{\epsilon _{x}}^{2}(X_{n}+k_{1x}\Delta t)^{2})\\k_{1y}=Y_{n}(-c+dX_{n}+{\epsilon _{y}}^{2}Y_{n}^{2})\\k_{2x}=(Y_{n}+k_{1y}\Delta t)(-c+d(X_{n}+k_{1x}\Delta t)+{\epsilon _{y}}^{2}(Y_{n}+k_{1y}\Delta t)^{2})\\\end{cases}}$

O índice $n$ representa o passo $\Delta t$ na iteração do método.

Stratonovich

O ruído também pode ser adicionado nos parâmetros $b$ e $d$, de forma que $b\to b+\epsilon _{x}\xi _{x}(t)$ e $d\to d+\epsilon _{y}\xi _{y}(t)$, onde o ruído segue a forma de Stratonovich. Essa mudança pode ser interpretada como a predação de diferentes predadores sobre uma mesma presa, considerando que não há competição entre as diferentes espécies. Esse ruído também pode ser pensado como uma variação do efeito da caça de um mesmo predador, ocasionado, por exemplo, a um certo número de presas escaparem à caça.

Dessa forma, a dinâmica populacional se comporta como:

${\begin{cases}{\dot {X}}=X(a-bY)-\epsilon _{x}XY\xi _{x}(t)\\{\dot {Y}}=Y(-c+fX)+\epsilon _{y}XY\xi _{y}(t)\end{cases}}$

onde $\xi (t)$ é um ruído branco e $\epsilon _{i}\geq 0$. Considerando que $\bullet$ denota a integral de Stratonovich e $W_{n}(t)$ um escalar oriundo de processos de Wiener independentes, utiliza-se para aplicação no método numérico a seguinte equação diferencial estocástica:

${\begin{cases}{\dot {X}}=X(a-bY)dt-\epsilon _{x}XY\bullet dW_{x}\\{\dot {Y}}=Y(-c+fX)dt+\epsilon _{y}XY\bullet dW_{y}\end{cases}}$

Utilizando o método de Heun tem-se:

${\begin{cases}X_{n+1}=X_{n}+{\frac {1}{2}}(A(x_{n})+A(X_{1}))\Delta t+{\frac {1}{2}}\Delta W_{x}(B(x_{n})+B(X_{1}))\\Y_{n+1}=Y_{n}+{\frac {1}{2}}(A(y_{n})+A(Y_{1}))\Delta t+{\frac {1}{2}}\Delta W_{y}(B(y_{n})+B(Y_{1}))\\\end{cases}}$

onde

${\begin{cases}A(x_{n})=x_{n}(a-by_{n})\\A(y_{n})=y_{n}(-c+fx_{n})\\\\B(x_{n})=\epsilon _{x}x_{n}y_{n}\\B(y_{n})=-\epsilon _{y}x_{n}y_{n}\\\\X_{1}=x_{n}+{\sqrt {\Delta t}}A(x_{n})+\Delta W_{n}B(x_{n})\\Y_{1}=y_{n}+{\sqrt {\Delta t}}A(y_{n})+\Delta W_{n}B(y_{n})\\\\A(X_{1})=X_{1}(a-bY_{1})\\A(Y_{1})=Y_{1}(-c+fX_{1})\\\\B(X_{1})=-\epsilon _{x}X_{1}Y_{1}\\B(Y_{1})=\epsilon _{y}X_{1}Y_{1}\\\end{cases}}$

Código

A estrutura do código foi a mesma para os três casos considerados, mudando apenas os detalhes da equação diferencial utilizada. A integração numérica da parte determinística foi feita pelo Método de Heun. Todas as análises foram realizadas a partir de um número N de realizações independentes, de onde se calcula a média e a variância para cada passo de tempo destas realizações. Assim, pode-se analisar as tendências do sistema de maneira estatística.

using LaTeXStrings
using Plots
using Random
using Statistics


module model
	function Lotka_Q_Ito(x,y,a,b,c,d,ϵ1,ϵ2,Δt) #Modelo Quadrático (Lucas)

		#Parte determinista da equação (Sem o ΔW)
		A_x = x*(a - b*y + (ϵ1^2)*x^2) 
		A_y = y*(-c + d*x + (ϵ2^2)*y^2)
		#Parte Estocástica
		B_x = ϵ1*x^2
		B_y = -ϵ2*y^2

		return A_x, A_y, B_x, B_y
	end
	function Lotka_Ito(x,y,a,b,c,d,ϵ1,ϵ2,Δt) #Modelo Linear Ito (Vini)
		#Parte determinista
		A_x = x*(a - b*y + (ϵ1^2)/2) 
		A_y = y*(-c + d*x + (ϵ2^2)/2)
		#Parte Estocástica
		B_x = ϵ1*x
		B_y = ϵ2*y
		return A_x, A_y, B_x, B_y
	end
	
	function Lotka_Strato_I(x,y,a,b,c,d,ϵ1,ϵ2,Δt) #Modelo Linear Stratonovich (André)
		#Parte determinista
		A_x = x*(a - b*y)
		A_y = y*(-c+d*x)
		#Parte Estocástica
		B_x = -ϵ1*x*y
		B_y = ϵ2*y*x
		return A_x, A_y, B_x, B_y
	end	
	
end

module simula

using LaTeXStrings
using Plots
using Random
using Statistics
using Main.model

	function Lotka(Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)

		n = floor(Int64, tmax/Δt) #Tamanho dos meus vetores
		rng = MersenneTwister() ##Chave para gerar números aleatórios
		ΔW_1 = sqrt(Δt)*randn(rng, n)
		ΔW_2 = sqrt(Δt)*randn(rng, n)
		
		x = zeros(n)
		y = zeros(n)
		x[1] = x0
		y[1] = y0
		
		i=1
		while i<n
			##Método de Heun (Runge Kutta 2)
			#Caso queira mudar o modelo utilizado, basta trocar o Lotka_Q_Ito pelo respectivo nas duas chamadas.

			k1_x, k1_y, B_x, B_y = model.Lotka_Ito(x[i],y[i],a,b,c,d,ϵ1,ϵ2, Δt)
			x_aux = x[i] + k1_x*Δt
			y_aux = y[i] + k1_y*Δt
			k2_x, k2_y = model.Lotka_Q_Ito(x_aux,y_aux,a,b,c,d,ϵ1,ϵ2, Δt)

			x[i+1] = x[i] +(k1_x + k2_x)*0.5*Δt + B_x*ΔW_1[i]
			y[i+1] = y[i] + (k1_y + k2_y)*0.5*Δt + B_y*ΔW_2[i]
			
			##Método de Heun (Stratonovich)
			#k1_x, k1_y, B_x, B_y = model.Lotka_Strato_I(x[i],y[i],a,b,c,d,ϵ1,ϵ2, Δt)
			#x_aux = x[i] + k1_x*sqrt(Δt) + B_x*ΔW_1[i]
			#y_aux = y[i] + k1_y*sqrt(Δt) + B_y*ΔW_2[i]
			#k2_x, k2_y, ka_x, ka_y = model.Lotka_Strato_I(x_aux,y_aux,a,b,c,d,ϵ1,ϵ2, Δt)

			#x[i+1] = x[i] +(k1_x + k2_x)*0.5*Δt + 0.5*(B_x+ka_x)*ΔW_1[i]
			#y[i+1] = y[i] + (k1_y + k2_y)*0.5*Δt + 0.5*(B_y+ka_y)*ΔW_2[i]
			
			if x[i+1]<0
				x[i+1]=0
			end
			if y[i+1]<0
				y[i+1]=0
			end

			##Método de Euler
			#x[i+1] = x[i] + k1_x*Δt + B_x*ΔW_1[i]
			#y[i+1] = y[i] + k1_y*Δt + B_y*ΔW_2[i]
			
			#Método de Euler Determinístico
			#x[i+1] = x[i] + x[i]*(a - b*y[i])*Δt 
			#y[i+1] = y[i] + y[i]*(-c + d*x[i])*Δt 
			
		i+=1
		end
		return x,y
	end
		
	function realiza(N, Δt,tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
		##Teste sobre N realizações
		n = floor(Int64, tmax/Δt)

		X = zeros(N, n) #Matrizes que armazenam valores de x para cada tempo em N realizações
		Y = zeros(N, n)
		mean_X = zeros(n)
		mean_Y = zeros(n)
		var_X = zeros(n)
		var_Y = zeros(n)


		i=1
		while i<=N
			X[i,:],Y[i,:] = simula.Lotka(Δt,tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
			i+=1
		end

		i=1
		while i<=n
			mean_X[i] = mean(X[:, i])
			mean_Y[i] = mean(Y[:, i])
			var_X[i] = var(X[:,i])
			var_Y[i] = var(Y[:,i])
			i+=1
		end

		return mean_X, mean_Y, var_X, var_Y
	end
end

module plotar

using LaTeXStrings
using Plots
using Main.simula


	function mapa(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
		##Essa função faz o mapa de fase para cada par de condições iniciais para x0 e y0
		##x0 e y0 devem ser vetores (mesmo que possuam somente um elemento), os outros todos devem ser escalares.
		##O plot é feito com a colorbar indicando a passagem de tempo conforme a progessão do tempo

		t = LinRange(0, tmax, floor(Int64, tmax/Δt))
		titulo = L"$a=$"*string(a) * L"$\: , b=$"*string(b) * L"$\: , c=$"*string(c) * L"$\:, d=$"*string(d) * L"$\:, \epsilon_1=$"*string(ϵ1) * L"$\: , \epsilon_2=$"*string(ϵ2)
		p=plot(title=titulo, xlabel="<X>", ylabel="<Y>", tickfontsize=11, guidefontsize=15, size=(1000,600))
		
		i=1
		while i<=length(y0)
			j=1
			while j<=length(x0)
				x, y, var_x, var_y = simula.realiza(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0[j], y0[i])
				legenda = L"$x_0=$"*string(x0[j])*L", $y_0=$"*string(y0[i])
				plot!(x, y, linez=t, color=:thermal,label=legenda, legendfont=font(15))
				j+=1
			end
			i+=1
		end
		display(p)
		savefig(p,"mapa_exemplo5.png")

	end

	function series(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
		##Plota a série temporal da média/variância de X sobre N realizações
		##Se quiser fazer para várias condições iniciais no mesmo gráfico fazer loop conforme mapa, mas a priori acho que não será necessário. Por isso, pego só o primeiro elemento de x0

		x, y, var_x, var_y = simula.realiza(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0[1], y0[1])

		n = floor(Int64, tmax/Δt)
		t = LinRange(0, tmax, n)

		titulo = L"$a=$"*string(a) * L"$\: , b=$"*string(b) * L"$\: , c=$"*string(c) * L"$\:, d=$"*string(d) * L"$\:, \epsilon_1=$"*string(ϵ1) * L"$\: , \epsilon_2=$"*string(ϵ2) *L"$\: ,x_0=$"*string(x0[1]) * L"$\: , y_0=$"*string(y0[1])
		
		px = plot(t, x, ylabel="<X>", xlabel=L"t", legend=false, color=:red)
		py = plot(t, y, ylabel="<Y>", xlabel=L"t", legend=false)
		px_var = plot(t, var_x, ylabel="Variância de X", xlabel=L"t", legend=false, color=:red)
		py_var = plot(t, var_y, ylabel="Variância de Y", xlabel=L"t", legend=false)
		
		p = plot(px, py, px_var, py_var, layout=(2,2), tickfontsize=11, guidefontsize=15, size=(1200, 600), plot_title=titulo)
		display(p)
		savefig(p,"series_exemplo5.png")

	end

	function anima(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
		t = LinRange(0, tmax, floor(Int64, tmax/Δt))
		titulo = L"$a=$"*string(a) * L"$\: , b=$"*string(b) * L"$\: , c=$"*string(c) * L"$\:, d=$"*string(d) * L"$\:, \epsilon_1=$"*string(ϵ1) * L"$\: , \epsilon_2=$"*string(ϵ2) *L"$\: ,x_0=$"*string(x0[1]) * L"$\: , y_0=$"*string(y0[1])
		p=plot(title=titulo, xlabel="<X>", ylabel="<Y>", tickfontsize=11, guidefontsize=15, size=(1000,600))
		

		x, y, var_x, var_y = simula.realiza(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0[1], y0[1])
		i=2
		anim = @animate while i<=length(t)
			#legenda = L"$x_0=$"*string(x0[j])*L", $y_0=$"*string(y0[i])
			titulo = L"$a=$"*string(a) * L"$\: , b=$"*string(b) * L"$\: , c=$"*string(c) * L"$\:, d=$"*string(d) * L"$\:, \epsilon_1=$"*string(ϵ1) * L"$\: , \epsilon_2=$"*string(ϵ2) *L"$\: ,x_0=$"*string(x0[1]) * L"$\: , y_0=$"*string(y0[1])
			p=plot(title=titulo, xlabel="<X>", ylabel="<Y>", tickfontsize=11, guidefontsize=15, xlim=(0,200), ylim=(0,300), size=(1000,600))
			plot!(x[1:i], y[1:i])
			#display(p)
			i+=1000
		end
		
		arquivo = "Nomedesejado.gif"
		gif(anim, arquivo, fps=500)

	end


end


function main() 

	##Condições Iniciais
	#Valores escolhidos conforme exemplo de paper
	a = 2.
	b = 0.5
	c = 0.5
	d = 1.

	Δt = 0.0001 ##Tem que ser pequeno mesmo
	tmax =200
	β = 1 #Fator de escala do ruído
	ϵ1 = β*0.001
	ϵ2 = β*0.001
	

	x0 = [10]#Número de presas iniciais	#equilíbrio: d/c
	y0 = [100]#Número de predadores iniciais	#equilibrio:a/b

	N = 50

	plotar.series(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0) #Se quiser plotar série basta trocar anima por serie, por exemplo.
	plotar.mapa(N, Δt, tmax, a, b, c, d, ϵ1, ϵ2, x0, y0)
end

main()

Uma outra versão do código, para o modelo linear de Itô, está disponível em https://github.com/Vini-BM/StochLV

Resultados

Nesta seção, serão abordados diversos aspectos sobre o comportamento numérico dos sistemas simulados. Avalia-se qualitativamente a estabilidade e o comportamento do sistema conforme a variação dos parâmetros $a,\ b,\ c,\ f$, das condições iniciais e também da intensidade do ruído.

Ruído externo

O modelo com ruído externo foi tratado apenas no sentido de Itô. As simulações foram feitas com o intuito de observar o efeito das condições iniciais e da intensidade do ruído.

Condições iniciais

Foram escolhidas algumas condições iniciais diferentes, com os mesmos parâmetros e ruído, a fim de observar o comportamento do sistema quanto ao ponto fixo. O espaço de fases e as séries temporais estão representados abaixo:

Espaço de fases para a condição inicial $(x^{\ast },y^{\ast })=(40.0,60.0)$

Espaço de fases para a condição inicial $(x^{\ast },y^{\ast })=(2.5,20.0)$

Séries temporais para a condição inicial $(x^{\ast },y^{\ast })=(40.0,60.0)$

Séries temporais para a condição inicial $(x^{\ast },y^{\ast })=(2.5,20.0)$

Desses resultados, percebe-se que o sistema estocástico não tem mais órbitas fechadas, e o ponto de equilíbrio se torna um atrator. Isso fica muito claro no primeiro mapa de fase, onde a trajetória do sistema até o ponto fixo é bastante direta e pouco ruidosa. Por outro lado, é interessante notar que o ponto de equilíbrio se torna instável, como visto nos gráficos da direita: ao começar em cima do ponto fixo, o ruído leva o sistema para uma "órbita" em torno dele. Observando as séries temporais, percebe-se que a média do sistema fica em torno do ponto de equilíbrio mas com uma variância alta, evidenciando o efeito do ruído sobre a estabilidade do sistema.

Intensidade do ruído

Passamos, então, ao estudo do ruído no sistema. Para tanto, fixamos os parâmetros e as condições iniciais e variamos a intensidade do ruído. Os resultados estão representados abaixo:

Espaço de fases para o ruído $\epsilon _{x}=0.01,\ \epsilon _{y}=0.02$

Espaço de fases para o ruído $\epsilon _{x}=0.1,\ \epsilon _{y}=0.2$

Séries temporais para o ruído $\epsilon _{x}=0.01,\ \epsilon _{y}=0.02$

Séries temporais para o ruído $\epsilon _{x}=0.1,\ \epsilon _{y}=0.2$

Observamos que o sistema se comporta relativamente bem com um ruído da ordem de $10^{-2}$, de acordo com as análises feitas anteriormente. Contudo, ao aumentar o ruído por um fator de $10^{1}$, a variância do sistema alcança valores maiores por um fator de $10^{3}$. O sistema ainda parece ter sua média em torno do ponto de equilíbrio, mas a variância é tão grande que o comportamento dele quase não é ditado por esse ponto.

Outro fator importante de se ressaltar é o fato de que $X$ e $Y$ não podem ser negativos. O ruído, sendo um processo de Wiener, é simétrico em torno do zero, portanto o sistema poderia divergir para $\pm \infty$. Porém, devido à condição $X,Y\geq 0$, que foi imposta sobre a integração numérica, o sistema dificilmente fica equilibrado em torno de zero, pois qualquer ruído negativo será descartado. Assim, para valores pequenos do ponto de equilíbrio, a média do sistema tem um comportamento assimétrico em torno dele. Quanto aos pontos de equilíbrio, isso significa que, enquanto o ponto fixo não trivial se torna uma espécie de ciclo limite, o ponto fixo trivial $(0,0)$ deixa de ser uma solução de equilíbrio.

Ruído nos parâmetros

O modelo com os parâmetros estocásticos foi analisado no sentido de Itô e no sentido de Stratonovich, conforme a dedução anterior.

Itô

Primeiramente realizou-se simulações variando o número de presas iniciais para um número fixo de predadores, representando o sistema em um mapa de fase da média de 50 realizações:

Nota-se um comportamento característico de atratores, onde o ponto fixo das equações de Lotka-Volterra determinística é o mesmo (ponto vermelho). Evidencia-se que esse ponto fixo não depende das condições iniciais.

O comportamento atrativo fica claro ao analisar a média e a variância em função do tempo:

As amplitudes para os predadores e as presas diminuem conforme o tempo, porém nota-se que o ruído passa a se tornar mais evidente quando as amplitudes das órbitas são menores.

Com a finalidade de testar a coerência do sistema com o determinístico, integrou-se as equações para um ruído da ordem de $10^{-9}$.

Mostrando, portanto, que o sistema estocástico se reduz ao determinístico no limite em que o ruído é muito pequeno, dado que se observam órbitas fechadas. Espera-se que a variância seja nula em um sistema determinístico, porém numericamente obteve-se oscilações que são da ordem do ruído ou menor, corroborando para a convergência do sistema estocástico ao determinístico com pequenos ruídos.

Ao realizar as simulações notou-se que o sistema é bastante sensível ao ruído, convergindo para valores bem estritos de parâmetros, ruído e condição inicial. Portanto, embora esse sistema modele uma interação entre dois níveis da cadeia alimentar, ele possui alguns graves problemas, como a fácil divergência -- devido, provavelmente, aos termos cúbicos das equações. Porém, também houve resultados interessantes como os atratores, que podem representar uma espécie de equilíbrio ecológico, onde não há extinções com o passar do tempo. Um dos problemas mais sérios com relação ao modelo é que mesmo na ausência de presas a população de predadores pode crescer indefinidamente, o que não condiz com a realidade.

As simulações realizadas variando os parâmetros trouxeram pequenas modificações no comportamento do sistema, como a alteração do período, a quantidade de órbitas bem definidas e suas amplitudes.

Nota-se que os períodos e as amplitudes diminuem conforme $a$ aumenta. Os ruídos se tornaram mais presentes com grandes variâncias nas simulações com $a$ pequeno. Outra tendência observada foi de que o regime ruidoso parece ocorrer mais tarde para valores de $a$ maiores, o que provavelmente está associado ao fato de o ruído se tornar mais presente para pequenas oscilações temporais do número de presas, o que ocorre na região próxima do ponto fixo.

Stratonovich

Com o objetivo de comparar a influência de diferentes ruídos com mesmos parâmetros, foi realizado o gráfico com condições iniciais iguais ao exemplo de Itô.

Percebe-se que as órbitas se distribuem de forma menos espaçada radialmente em relação à simulação de Itô e possuem trajetórias menos ruidosas.

Em seguida, fez-se o gráfico do comportamento da média e variância das densidades populacionais. Pode-se observar que os valores médios convergem para os valores do ponto fixo e continuam variando em torno dele numa forma de batimento. A amplitude das variações da variância não se estabilizou.

No próximo teste, variou-se o valor do ruído procurando um valor máximo para a simulação não divergir. Encontrou-se o valor $\epsilon =0.00046$, porém os sistema pode não convergir para tempos de integração muito grandes. Observou-se que isso independe da condição inicial.

Por fim, a variação dos parâmetros $a$, $b$, $c$ e $d$ pode comprometer a convergência da simulação, como, por exemplo, para os casos de taxa de crescimento do predador muito alta. Simulou-se as equações com parâmetros distintos para evidenciar o comportamento bastante convergente do método em condições iniciais extremas:

Referências

1. BRAUER, F.; CASTILLO-CHAVEZ, C. Mathematical Models in Population Biology and Epidemiology. New York, NY: Springer New York, 2012. v. 40
2. COELHO, P. J. de O. Equações de Lotka–Volterra Estocásticas: Simulações com o Matlab. 25 de junho de 2015. Disponível em https://www.academia.edu/52185574/Sistema_de_competi%C3%A7%C3%A3o_Lotka_Volterra_sob_ru%C3%ADdo_branco. Acesso em ago. 2024.
3. SCHERER, C. Métodos computacionais da física. 2. ed ed. São Paulo: Liv. da Física, 2010.
4. KHASMINSKII, R. Z.; KLEBANER, F. C. Long term behavior of solutions of the Lotka-Volterra system under small random perturbations. The Annals of Applied Probability, v. 11, n. 3, 1 ago. 2001.
5. KOERS, L. Geometric integration of stochastic Lotka-Volterra equations. BSc Thesis Applied Mathematics, 8 ago. 2024.