Numerical Methods for Solving the Lotka–Volterra Equations

Introduction

The Lotka–Volterra system arises in mathematical biology and models the growth of animal species. Consider two species where (Subscript[Y, 1])(T) denotes the number of predators and (Subscript[Y, 2])(T) denotes the number of prey. A particular case of the Lotka–Volterra differential system is

({{(OverDot[(Subscript[Y, 1])]), ==, (Subscript[Y, 1])((Subscript[Y, 2])-1),, (OverDot[(Subscript[Y, 2])]), ==, (Subscript[Y, 2])(2-(Subscript[Y, 1]))}}),

where the dot denotes differentiation with respect to time T.

The Lotka–Volterra system () has an invariant H, which is constant for all T:

H((Subscript[Y, 1]),(Subscript[Y, 2]))==2ln (Subscript[Y, 1])-(Subscript[Y, 1])+ln (Subscript[Y, 2])-(Subscript[Y, 2]).

The level curves of the invariant () are closed so that the solution is periodic. It is desirable that the numerical solution of () is also periodic, but this is not always the case.  Note that () is a Poisson system:

({{(OverDot[Y]), ==, B(Y) ∇H(Y), ==, (({{0, -(Subscript[Y, 1])(Subscript[Y, 2])}, {(Subscript[Y, 1])(Subscript[Y, 2]), 0}}))(({{(2/(Subscript[Y, 1]))-1}, {(1/(Subscript[Y, 2]))-1}}))}})

where H(Y) is defined in ().

Poisson systems and Poisson integrators are discussed in Chapter VII.2 of [HLW02] and [MQ02].

Load a package with some predefined problems and select the Lotka–Volterra system:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

system = GetNDSolveProblem["LotkaVolterra"];
invts = system["Invariants"];
time = system["TimeData"];
vars = system["DependentVariables"];
step = 3/25;

Define a utility function for visualizing solutions:

LotkaVolterraPlot[sol_,vars_,time_,opts___?OptionQ]:=
Module[{commonopts, data,data1,data2,ifuns,lplot,pplot},
ifuns=First[vars/.sol];
data1=Part[ifuns,1,0]["ValuesOnGrid"];
data2=Part[ifuns,2,0]["ValuesOnGrid"];
data=Transpose[{data1,data2}];
commonopts=Sequence[Axes->False,Frame->True,FrameLabel->Join[Map[TraditionalForm,vars],{None,None}],RotateLabel->False];
lplot=ListPlot[data,Evaluate[FilterRules[{opts},Options[ListPlot]]],PlotStyle->{PointSize[0.02],RGBColor[0,1,0]},Evaluate[commonopts]];
pplot=ParametricPlot[Evaluate[ifuns],time,Evaluate[FilterRules[{opts},Options[ParametricPlot]]],Evaluate[commonopts]];
Show[lplot,pplot]
];

Explicit Euler

Use the explicit or forward Euler method to solve the system ():

fesol = NDSolve[system,Method->"ExplicitEuler",StartingStepSize->step];

LotkaVolterraPlot[fesol, vars, time]

Backward Euler

Define the backward or implicit Euler method in terms of the RadauIIA implicit Runge–Kutta method and use it to solve (). The resulting trajectory spirals from the initial conditions toward a fixed point at (2,1) in a clockwise direction:

BackwardEuler ={"FixedStep",Method-> {"ImplicitRungeKutta","Coefficients"->"ImplicitRungeKuttaRadauIIACoefficients","DifferenceOrder"->1,"ImplicitSolver"->{"FixedPoint",AccuracyGoal->MachinePrecision, PrecisionGoal->MachinePrecision,"IterationSafetyFactor"->1}}};

besol = NDSolve[system,Method->BackwardEuler,StartingStepSize->step];

LotkaVolterraPlot[besol, vars, time]

Projection

Projection of the forward Euler method using the invariant () of the Lotka–Volterra equations gives a periodic solution:

pfesol = NDSolve[system,Method->{Projection,Method->"ExplicitEuler", Invariants->invts},StartingStepSize->step];

LotkaVolterraPlot[pfesol, vars, time]

Splitting

Another approach for obtaining the correct qualitative behavior is to additively split () into two systems:

({{(OverDot[(Subscript[Y, 1])]), ==, (Subscript[Y, 1]) ((Subscript[Y, 2])-1 ), (OverDot[(Subscript[Y, 2])]), ==, 0}, {(OverDot[(Subscript[Y, 1])]), ==, 0, (OverDot[(Subscript[Y, 2])]), ==, (Subscript[Y, 2]) (2-(Subscript[Y, 1])).}})

By appropriately solving () it is possible to construct Poisson integrators.

Define the equations for splitting of the Lotka–Volterra equations:

eqs = system["System"];
Y1 = eqs;
Part[Y1, 2, 2] = 0;
Y2 = eqs;
Part[Y2,1, 2] = 0;

Symplectic Euler

Define the symplectic Euler method in terms of a splitting method using the backward and forward Euler methods for each system in ():

SymplecticEuler= {"Splitting",
 "DifferenceOrder"->1,"Equations"->{Y1,Y2}, "Method"->{BackwardEuler,"ExplicitEuler"}};

sesol = NDSolve[system,Method->SymplecticEuler,StartingStepSize->step];

The numerical solution using the symplectic Euler method is periodic:

LotkaVolterraPlot[sesol, vars, time]

Flows

Consider splitting the Lotka–Volterra equations and computing the flow (or exact solution) of each system in (). The solutions can be found as follows, where the constants should be related to the initial conditions at each step.

DSolve[Y1,vars,T]

DSolve[Y2,vars,T]

An advantage of locally computing the flow is that it yields an explicit, and hence very efficient, integration procedure. The "LocallyExact" method provides a general way of computing the flow of each splitting using (RefLink[DSolve, paclet:ref/DSolve]) only during the initialization phase.

Set up a hybrid symbolic-numeric splitting method and use it to solve the Lotka–Volterra system:

SplittingLotkaVolterra= {"Splitting",
 "DifferenceOrder"->1,"Equations"->{Y1,Y2}, "Method"->{"LocallyExact","LocallyExact"}};

spsol = NDSolve[system,Method->SplittingLotkaVolterra,StartingStepSize->step];

The numerical solution using the splitting method is periodic:

LotkaVolterraPlot[spsol, vars, time]