170 likes | 278 Views
A Simple Model of Carnivore & Their Prey. How to deal with linked ODE systems And why we must abandon Euler's method Next time; how to replace Euler's. First, what is an ODE? Where the derivatives of one or more variables vary with respect to just one axis
E N D
A Simple Model of Carnivore & Their Prey How to deal with linked ODE systems And why we must abandon Euler's method Next time; how to replace Euler's
First, what is an ODE? • Where the derivatives of one or more variables vary with respect to just one axis • So let us examine a simple model of how Carnivores and Prey evolve in time (but not space & time). • Why is this earth science?
Lotka-Volterra Predator Prey model • Assume prey have infinite food, and each can reproduce a fixed amount per time • So prey (P) growth is a fraction of prey abundance • This really means that prey abundance is controlled by carnivores (why?) • Solution, without predators, is • Show on board!
Assume that a fraction of the carnivores C die per unit time • accidents and mishaps! • so w/o food, they all die in an exponential decay
And now the interaction term • number of Prey eaten depends on number of prey times number of carnivores • Chance encounter model • with efficiency • thus: • Where PC is rate of predation
And some of the Prey eaten turn into new carnivores with an efficiency • burp…
So the equations are: Prey are born Prey are eaten, and turned into Carnivores Carnivores die of old age
How do we solve this with Euler's method? • So we write code for the derivative of the vector Pop=[P,C]
So what does the code for the derivative function look like? function dPopdt=dPop_dt(t,Pop); %calculate the Lotka-Volterra carnivore prey equations % % d(Prey)/dt = Prey*(alpha-beta*Carnivore) % d(Carnivore)/dt = -Carnivore*(gamma-delta*Prey) % % where % % The state variable is Pop(1) is prey, Pop(2) is carnivore gamma=0.05; delta=0.001; alpha=0.04; beta=0.00075; dPopdt=[1;1]; %Explain this line! dPopdt(1)=Pop(1)*(alpha-beta*Pop(2)); dPopdt(2)=-Pop(2)*(gamma-delta*Pop(1));
And how do we use this with Euler's method? t=0; %initial time Pop=[10;10]; %intial condition dt=0.05/2; %time step tvec=[t]; %store time history Popvec=[Pop]; %store T history %integrate with Euler's method to time t>2pi while (t<=500) Pop=Pop+dt*dPop_dt(t,Pop); %what is the shape of Pop? t=t+dt; tvec=[tvec t]; Popvec=[Popvec Pop]; end plot(tvec,Popvec(2,:),'r*-',tvec,Popvec(1,:),'g-','linewidth',2) xlabel('time, days') ylabel('critter numbers') title('Red line is the Carnivore density; green is the Prey')
So what does the solution look like • What is wrong? • What is right?
What is right: The dynamics • Few Carnivores • Lots of Prey growth • Lots of Preditor growth • Most of the prey eaten • Carnivores Starve
Examples of this in the real world • diatom blooms in ocean • lemmings in the arctic • wolf & moose in some isolated systems (islands, mostly)
What is wrong? • With Euler, we do not know what timestep should be
How do you find a good time step? • What is origin of error in Euler's?
So how do you find the "right" timestep? • should get same answer for t and t/2 • But this is expensive and time consuming! • And Euler's is also inefficient. • Next time, a more accurate method that chooses its own timestep!