Gradacją w leśnictwie nazywa się wystąpienie dużej liczby owadów w danym sezonie. Celem modelu jest zbadanie warunków pojawienia się gradacji. Model rozwoju szkodnika drzew iglastych (głównie jodły balsamicznej) z gatunku Choristoneura fumiferana (Clemens) – larwy motyla z rodziny zwójkowatych zaproponowany w 1978 roku przez D. Ludwiga. Równanie opisujące populację szkodników wygląda następująco:
gdzie:
N – populacja szkodników
rB – współczynnik rozrodczości
KB – pojemność środowiska związana z ilością igieł na drzewach, czyli pożywienia dla owadów
p(N) – funkcja opisująca drapieżnictwo ze strony ptaków, dla których owady te są pożywieniem, funkcją p(N) definiujemy następująco:
parametry A i B określają drapieżnictwo.
Typowy przedstawiciel gatunku Choristoneura fumiferana wygląda tak na na ilustracji poniżej:
Ilustracja
1: Choristoneura fumiferana. Autor: Jerald
E. Dewey
http://www.forestryimages.org/browse/detail.cfm?imgnum=2252020
Model wykonany w simulinku wygląda następująco:
Ilustracja
2: model w simulinku
zaś funkcja p(N):
Ilustracja
3: funkcja p(N)
Model pozwala nam na narysowanie populacji owadów N w czasie, oraz zobaczyć postać funkcji p(N). W zależności od wartości parametrów opisujących model możemy otrzymać wzrost populacji, jej zmniejszenie, lub brak zmian.
Seria ilustracji poniżej przedstawia przykładowe wyniki:
Ilustracja
4: spadek liczebności populacji
Ilustracja
5: funkcja p(N) odpowiadająca populacji z ilustracji 4
Ilustracja
6: wzrost liczebności populacji
Ilustracja
7: funkcja p(N) odpowiadająca populacji z ilustracji 6
Ilustracja
8: brak znaczący zmian w liczbie populacji
Ilustracja
9: funkcja p(N) odpowiadająca populacji z ilustracji 8
Ten sam problem można rozwiązać w Matlabie.
function F = modelGradacji(t,N)
% N – populacja szkodnikow
% t - czas
rB = 0.6; % wspolczynnik rozrodczosci
KB = 0.5; % pojemnosc srodowiska zwiazana z iloscia igiel na drzewach
% funkcja opisujaca drapieznictwo ze strony ptakow
A = sqrt(0.3);
B = 4;% parametry A i B okreslaja drapieznictwo.
F = rB*N.*(1-N/KB)-B*N.^2./(A^2+N.^2);
Przy użyciu wbudowanego solvera, np. ode45:
[t,y]= ode45('modelGradacji',[0 10],0.1);
plot(t,y)
xlabel('czas')
ylabel('liczebnosc N')
title('populacji szkodnikow')
Ilustracja
10: rozwiązanie rownania za pomocą metody ode45
lub za pomocą 'ulepszonej' metody Eulera:
function [t,y]=euler_ulepszona(f,tinit,yinit,tfinal,n)
% Ulepszona metoda Eulera
% Obliczenia kroku
h=(tfinal-tinit)/n;
% Przygotowanie wektorów poczatkowych t i y
t=[tinit; zeros(n,1)]; y=[yinit; zeros(n,1)];
% Obliczenia t i y
for i=1:n
t(i+1)=t(i)+h;
ynew=y(i)+h*f(t(i),y(i));
y(i+1)=y(i)+(h/2)*(f(t(i),y(i))+f(t(i+1),ynew));
end
[t y]=euler_ulepszona(@modelGradacji,0,0.1,10,100);
figure
plot(t,y)
xlabel('czas')
ylabel('liczebnosc N')
title('populacja szkodnikow metoda: euler ulepszona')
Ilustracja
11: rozwiazanie za pomoca ulepszonej metody Eulera
Źródło: http://home.agh.edu.pl/~zygmunt/Modele/wyklad_modele_3.pdf
Pliki do pobrania: gradacja.zip
Grzegorz Knor 2010