function test9 x = [ 29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.3 69.6 39.8 34.0 ... 20.7 21.7 37.6 57.6 124.6 215.8 272.7 195.7 95.0 41.9 25.7 ... 10.9 22.6 33.6 48.1 92.5 183.3 268.5 230.6 111.1 ]'; y = [ 128 104 88 96 88 104 144 176 192 168 152 ... 120 120 96 72 88 104 128 184 192 168 136 ... 144 112 96 80 72 88 104 152 184 ]'; n = size(x,1); plot(1:n,x','r-',1:n,y','b-','linewidth',2); legend('mice','owl'); axis([0,50,0,300]); pause(3) yb = ( y(1:n-1)+y(2:n) )/2; xb = ( x(1:n-1)+x(2:n) )/2; % yb = [ (y(1)+y(2))/2; (y(1:n-3)+3*y(2:n-2)+3*y(3:n-1)+y(4:n))/8; (y(n-1)+y(n))/2 ]; % xb = [ (x(1)+x(2))/2; (x(1:n-3)+3*x(2:n-2)+3*x(3:n-1)+x(4:n))/8; (x(n-1)+x(n))/2 ]; A = [ +ones(n-1,1) -yb ]; B = [ -ones(n-1,1) +xb ]; a = diff(log(x)); b = diff(log(y)); p = A \ a; q = B \ b; lambda = p(1); alpha = p(2); mu = q(1); beta = q(2); [t,Y] = ode45(@lotka,[0 50],[x(1) y(1)]',odeset,lambda,mu,alpha,beta); close; figure('position',[50 50 1200 600]) subplot(1,2,1); plot(1:n,x','r-',t,Y(:,1),'k-','linewidth',2); legend('real','computed'); axis([0,50,0,300]); title('mice','fontsize',15); subplot(1,2,2); plot(1:n,y','r-',t,Y(:,2),'k-','linewidth',2); legend('real','computed'); axis([0,50,0,300]); title('owl','fontsize',15); pause(3); close; plot(Y(:,1),Y(:,2),'r-','linewidth',2); function dydx = lotka(t,y,lambda,mu,alpha,beta) dydx = [ lambda*y(1)-alpha*y(1)*y(2) -mu*y(2)+beta*y(1)*y(2) ];