Can someone explain to me whats off with this function handle?

So I have a function where I want to model the population of predators and prey based off of the classic model Lotka-Volterra. In my code I designate dy as an array with the two equations.

function[predator,prey]=lotka_volterra(fprime, timespan, y0, h)

h=.1; %step size
%initial conditions
X(1)=timespan(1);
Y(1)=4;
Y(2)=4;
prey=Y(1);
predator=Y(2);
dy=[2*prey(X)-predator(X)*prey(X); predator(X)*prey(X)-2*predator(X)];

i=2;
while X(end)<timespan(end)
    X(i)=X(i-1) + h;
    Y(i)=Y(i-1)+h*fprime(X(i-1)); % Y(i)=Y(i+1)+h*y'(X(i-1))
    i=i+1;

end

   end

In the command window I typed,

[predator,prey]=lotka_volterra(@(X)(dy), [0 10], [4 4], .1)

but it is not recognizing that the equations are supposed to change by X (time). What am I messing up?

[Andrew Newell]

It probably has something to do with the fact that you don't use dy after you define it.

[Jan]

"Is not recognized" is a vague description of the problem only. Do you get an error message? Then please read it carefully and if this does not reveal the details for you, please post a complete copy of the message in the forum.

[Adam]

What are you expecting

@(X)(dy)

to do? dy would have to be defined in the workspace above that call for this to make even some sense, but the placeholder 'X' is not used at all so has no purpose.

[Mohannad Abboushi]

So I changed the code a bit:

function[X,Y]=lotka_volterra_euler(yprime,timespan, y0, h)

h=.1; %step size
%initial conditions
X=[];
X(end+1)=0;
X(end+1)=timespan(end);
Y(1)=4;
Y(2)=4;

dx=2*Y(1)-Y(1)*Y(2);
dy=2*Y(1)*Y(2)-2*Y(2);
yprime=[dx dy];
yprime = @(X, Y) [dx dy];

if size(Y,1)>1; Y=Y'; end %if y0 is not scalar, make sure to keep it as a row vector.

i=2;
while true
    if (X(i-1)+h) > timespan(end); break; end
    nextX=X(i-1)+h;
    X(i)=X(i-1) + h;
  nextY = Y(i-1,:) + h*yprime(X(i-1), Y(i-1,:)); %nextY can be a scalar or a vector value.
    X(i) = nextX;
  Y(i,:) = nextY;
  i=i+1;
end

figure
plot(X,Y(:,1),X,Y(:,2));
xlabel('Years');
ylabel('Population (in thousands)');
title('Lotka-Volterra Graph');
legend('Euler approximation');

end 
Now when I put into the command window: [X,Y]=lotka_volterra_euler( @(X, Y) [dx dy],[0 10], [4 4], .1) I get a figure but the Y columns are not behaving cyclically as a population model would.

Answer 1

[Steven Lord]

This code doesn't do what you think it does.

dx=2*Y(1)-Y(1)*Y(2);
dy=2*Y(1)*Y(2)-2*Y(2);
yprime=[dx dy];
yprime = @(X, Y) [dx dy];

It computes NUMERIC values for dx and dy based on the values currently stored in the variable Y. When you define yprime the first time, it concatenates those numeric values together. When you define yprime the second time, it creates a function handle that will return those numeric values. Note that even though the function handle yprime accepts X and Y as inputs, dx and dy are numbers, not functions of X and Y, and so will not change. This function handle will return the same values every time, ignoring the values with which you call it.

You can do this using multiple function handles:

dx = @(X, Y) 2*X-X.*Y;
dy = @(X, Y) 2*X.*Y-2*Y;
yprime = @(X, Y) [dx(X, Y), dy(X, Y)];

Now when you evaluate yprime, it evaluates the function handles dx and dy and uses the output from those two function handles to generate the output that yprime will return.

[Mohannad Abboushi]

Awesome that makes sense! However when I use that formatting, I get the error:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in lotka_volterra_euler (line 23)
  Y(i,:) = nextY;

Answer 2

[Jan]

Try this:

[predator, prey] = lotka_volterra(@dy, [0 10], [4 4], 0.1)