MODELOWANIE BIOSYSTEMÓW

Laboratorium 2.

Temat ćwiczenia: Modelowanie pojedynczych populacji.

Zadanie 1.

a) kod programu zoptymalizowany w porównaniu z napisanym w trakcie laboratorium, wygładzone wykresy dla modelu ciągłego.

close all;

clear all;

r1 = [0.3, 1, 1.5]; % wektor wartoœci r_1

r2 = [0.3, 1, 1.5]; % wektor wartoœci r_2

t_d = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; % wektor czasu dla modelu dyskretnego

t_c = 0:.1:9; % wektor czasu dla modelu ciaglego

% model ciagly dla N_0 = 5

for i = 1:3

N_c1(i, 1) = 5;

for j = 1:90

N_c1(i, j+1) = N_c1(i, 1)*exp(r1(i)*t_c(j));

end;

end;

% model ciagly dla N_0 = 8

for i = 1:3

N_c2(i, 1) = 8;

for j = 1:90

N_c2(i, j+1) = N_c2(i, 1)*exp(r1(i)*t_c(j));

end;

end;

% model dyskretny dla N_0 = 5

for i = 1:3

N_d1(i, 1) = 5;

for j = 1:9

N_d1(i, j+1) = N_d1(i, j)*(1+r2(i));

end;

end;

% model dyskretny dla N_0 = 8

for i = 1:3

N_d2(i, 1) = 8;

for j = 1:9

N_d2(i, j+1) = N_d2(i, j)*(1+r2(i));

end;

end;

% wykresy

for i = 1:3

figure(i);

plot(t_d, N_d1(i, :), 'o', t_d, N_d2(i, :), 'ro', t_c, N_c1(i, :), 'g-', t_c, N_c2(i, :), 'k-');

legend('dyskretny dla N_0 = 5', 'dyskretny dla N_0 = 8', 'ciagly dla N_0 = 5', 'ciagly dla N_0 = 8');

xlabel('os czasu'); ylabel('N');

end;

b) wartość r₁ podana przez prowadzącego: r₁ = 0,3

dN(t)/dt = r₁ * N(t)

N(t) = N₀ * exp(r₁ * t)

N(t = T_D) = 2*N₀

N(T_D) = N₀ * exp(r₁ * T_D) = 2*N₀

exp(r₁ * T_D) = 2

r₁ * T_D = ln 2

T_D = ln 2/ r₁

T_D = 2,31

100. wartości podane przez prowadzącego: N₀ = 3, r₁ = 1,2;

N(t) = N₀ * exp(r₁ * t)

N(t)/ N₀ = exp(r₁ * t)

r₁ * t = ln(N(t)/ N₀) /:t

r₁ =

r₁ = 0,78845

500. wartość r₂ podana przez prowadzącego: r₂ = 1,2

Dla modelu dyskretnego: N_t = N_t-1 * (1+r₂) = N₀ * (1+r₂)^t

Dla modelu ciągłego: N(t) = N₀ * exp(r₁ * t)

N₀ * exp(r₁ * t) = N₀ * (1+r₂)^t /:N₀

e^{r1 * t} = (1+r₂)^t / ln

r₁ * t = t * ln(1+r₂) /:t

r₁ = ln(1+r₂)

r₁ = 0,78845

clear all;

close all;

r_2 = 1.2; % podane r_2

r_1 = log(1+r_2); % wyliczone r_1

t_d = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; % wektor czasu dla modelu dyskretnego

t_c = 0:.1:9; % wektor czasu dla modelu ciągłego

N_0 = 3;

% model ciagly

N_c = N_0.*exp(r_1*t_c);

% model dyskretny

N_d = N_0.*((1+r_2).^t_d);

% wykres

plot(t_c, N_c, t_d, N_d, 'ro');

legend('model ciagly', 'model dyskretny');

xlabel('os czasu'); ylabel('N');

Zadanie 2.

function derivate_N = logistyczny1(t, N)

a = 0.1;

N_infinity = 500;

derivate_N = a.*N.*(1-(N./N_infinity));

end

function derivate_N = logistyczny2(t, N)

a = 0.2;

N_infinity = 500;

derivate_N = a.*N.*(1-(N./N_infinity));

end

close all;

clear all;

zakres_t = [0, 120];

N_0 = [5, 50, 400, 750];

[t, N] = ode45(@logistyczny1, zakres_t, N_0);

plot(t, N(:,1), 'r', t, N(:,2), 'g', t,N(:,3),'b', t,N(:,4),'k');

xlabel('czas'); ylabel('ilosc komorek');

title('Przebieg czasowy modelu logistycznego dla a = 0,1');

legend('N_0 = 5','N_0 = 50','N_0 = 400','N_0 = 750');

grid on;

function deriv_N = gompertz(t, N)

b = 0.1;

N_infinity = 500;

deriv_N = -b.*N.*log(N./N_infinity);

end

zakres_t = [0, 120];

N_0 = [5, 50, 400, 750];

[t, N] = ode45(@gompertz, zakres_t, N_0);

plot(t, N(:,1), 'r', t, N(:,2), 'g', t,N(:,3),'b', t,N(:,4),'k');

xlabel('czas'); ylabel('ilosc komorek');

title('Przebieg czasowy modelu Gompertza dla b = 0,1');

legend('N_0 = 5','N_0 = 50','N_0 = 400','N_0 = 750'); grid on;

close all;

clear all;

zakres_t = [0, 120];

N_0 = [5, 50, 400, 750];

figure(1);

for i = 1:4

[T, N] = ode45(@logistyczny1,zakres_t ,N_0(i));

title('Wykresy fazowe dla modelu logistycznego, a = 0,1');

subplot(2,2,i);

plot(N(1:length(N)-1), N(2:length(N)));

end;

figure(2);

for i = 1:4

[T, N] = ode45(@logistyczny2,zakres_t ,N_0(i));

title('Wykresy fazowe dla modelu logistycznego, a = 0,2');

subplot(2,2,i);

plot(N(1:length(N)-1), N(2:length(N)));

end;

figure(3);

for i=1:4

[T, N] = ode45(@gompertz, zakres_t, N_0(i));

title('Wykresy fazowe dla modelu Gompertza, b = 0,1');

subplot(2,2,i);

plot(N(1:length(N)-1),N(2:length(N)));

end;

Zadanie 3.

a)

close all;

clear all;

a = [1.3, 2.5, 3.2, 3.5, 3.7]; % wektor wartosci a

t = [0:60]; % wektor czasu

X(1:5, 1) = 0.1; % podana wartosc X_0

for i = 1:5 % zmieniamy a

for j = 1:60

X(i, j+1)=a(i).*X(i, j).*(1-X(i, j));

end;

end;

for i = 1:5

figure(i);

plot(t, X(i, :), '*', t, X(i, :), 'g-');

xlabel('czas'); ylabel('X(t)');

title('Przebieg czasowy modelu logistycznego dyskretnego');

end;

dla a = 1,3

dla a = 2,5

dla a = 3,2

dla a = 3,5

dla a = 3,7

b) Oscylacje występują dla a = 3,2. Chaos dla a = 3.7.

c) diagram bifurkacyjny

bifurkacyjny.m

N = 50;

N1 = 25;

s = zeros(1, N);

s(1) = 0.1;

for j = 2:N

s(j) = a*s(j-1)*(1-s(j-1));

end;

figure(1);

plot(s,'o-');

figure(2);

plot(s(N1:end-1),s(N1+1:end),'.');

bif = [];

for a = 0:0.01:4

bifurkacyjny;

drawnow;

bif(end+1,:) = s(N1:end);

end;

figure();

plot([0:0.01:4],bif,'.','color','g','markersize',2);

---