How to plot DE system if functions are "periodic"?

Question

I use the Lotka–Volterra predator–prey model, so I have a system of two DE, and their solution are "periodic" functions. My code is given below.

with(DEtools): 
de_sysset := diff(x(t), t) = -.1*x(t)+0.2e-1*x(t)*y(t), diff(y(t), t) = .2*y(t)-0.25e-1*x(t)*y(t); 
de_sys := [de_sysset]; 
DEplot(de_sys, [x(t), y(t)], t = 0 .. 100, {[x(0) = 6, y(0) = 6]});

DEplot is correct, so I want to plot x(t) and y(t) on the same plot. I've started with

sol_sys := dsolve({de_sysset, x(0) = 6, y(0) = 6}, {x(t), y(t)}, numeric):

but I don't know how to work with plots from this point. All my attempts to plot give different errors. Using (Matlab) Simulink gives right result, so I think that everything is OK with my DE system. How can I plot solution of system of DE? Also, how can I plot all solutions together if they are for diffirent initial conditions (say, I have second condition [x(0) = 32, y(0) = 24]). Besides, I'm interested if I can obtain implicit solutions for DEs somehow.

I've also noticed that for large t intervals (say, t=4000), DEplot becomes "meshy", and it's more "meshy" in its upper right corner and more smooth in lower left. I wonder what causes that.

Answer 1

You can do these sort of plots using DEtools:-DEplot (passed the DE), or using plots:-odeplot (passed the result from dsolve,numeric), or using individual procedures (extracted from the result from dsolve,numeric with output=listprocedure).

See the help pages for more information on the options that can be used with these commands.

First, your original,

restart;
with(DEtools):
de_sysset := diff(x(t), t) = -.1*x(t)+0.2e-1*x(t)*y(t), diff(y(t), t) = .2*y(t)-0.25e-1*x(t)*y(t):
de_sys := [de_sysset]:

DEplot(de_sys, [x(t), y(t)], t = 0 .. 100, {[x(0) = 6, y(0) = 6]});

I prefer working with output=listprocedure, and using 2-argument eval to get at the solutions procedures for the dependent variables individually.

sol_sys := dsolve({de_sysset, x(0) = 6, y(0) = 6}, {x(t), y(t)}, numeric, output=listprocedure);

        sol_sys := [t = proc(t)  ...  end;, x(t) = proc(t)  ...  end;, 
                     y(t) = proc(t)  ...  end;]

X := eval(x(t), sol_sys);

                X := proc(t)  ...  end;

Y := eval(y(t), sol_sys);

                Y := proc(t)  ...  end;

# parametric form, agreeing with earlier DEplot
plot([X(t), Y(t), t=0..100]);

# separately as functions of t
plot([X(t), Y(t)], t=0..100, color=["Burgundy", "Navy"]);

Now, using plots:-odeplot,

plots:-odeplot(sol_sys, [x(t), y(t)], t=0..100);

plots:-display(
  plots:-odeplot(sol_sys, [t, x(t)], t=0..100, color="Burgundy"),
  plots:-odeplot(sol_sys, [t, y(t)], t=0..100, color="Navy")
);

You can also use DEtools[DEplot] to get the individual curves for x(t) or y(t). This may not be as efficient, and the plotting is not adapative so a larger numpoints value can help make the curve be smoother.

plots:-display(
  DEplot([de_sysset], [x(t), y(t)], t = 0 .. 100,
         scene=[t, x(t)], {[x(0) = 6, y(0) = 6]},
         linecolor="Niagara Burgundy",
         numpoints=100, thickness=1),
  DEplot([de_sysset], [x(t), y(t)], t = 0 .. 100,
         scene=[t, y(t)], {[x(0) = 6, y(0) = 6]},
         linecolor="Niagara Navy",
         numpoints=100, thickness=1)
);

Now, for multiple pairs of initial conditions, with some custom colors,

inits := { [x(0) = 6,  y(0) = 6],   [x(0) = 14, y(0) = 14],
           [x(0) = 20, y(0) = 18],  [x(0) = 26, y(0) = 22],
           [x(0) = 32, y(0) = 24] }:

colorlist := [seq(ColorTools:-Color([0,i,1-i]),
                  i=0.1 .. 0.9, 1/nops(inits))]:

DEplot(de_sys, [x(t), y(t)], t = 0 .. 100, inits,
       numpoints=1000, thickness=2, linecolor=colorlist);

Note that DEplot doesn't automatically adjust the number of computed points which make up the curve. So if you make the end-time very large then you get a more jagged curve from the default number of computed points.

DEplot(de_sys, [x(t), y(t)], t = 0 .. 4000, {[x(0) = 6, y(0) = 6]},
       thickness=0, linecolor=blue);

By increasing the number of computed points we can get back to a smoother curve, for the larger end-point.

DEplot(de_sys, [x(t), y(t)], t = 0 .. 4000, {[x(0) = 6, y(0) = 6]},
       numpoints=4000, thickness=0, linecolor=blue);