High-precision Approximate Closed Form Solutions For
 Nonlinear 2×2 Systems Of First-order IVPs  
 With Periodic Solutions

Temple H Fay, Technikon Pretoria and the University of Southern Mississippi
Stephan V Joubert and Alun J Peck, Technikon Pretoria, South Africa
joubertsv@hotmail.com

Introduction

This notebook provides an interactive routine for producing high-precision approximate solutions (truncated Fourier Series) for certain low-degree nonlinear 2×2  dynamical systems of first order ODEs (IVPs) with periodic solutions.  The routine is primarily meant for use in the classroom by undergraduate and graduate students, but might also be useful to technologists and scientists working with such systems for which no obvious closed-form solution exists.  It has the advantage that only the coefficients and the period need be stored and then the solution can be transferred to another platform such as a hand held calculator.

The 2×2 systems we will be considering are of the form:

x'[t] - f[x, y, t] = 0
y'[t] - g[x, y, t] = 0
x[0]=α   and    y[0]=β

and which are known to have periodic solutions . Furthermore, each system considered must be of such a nature that the chosen NDSolve routine performs satisfactorily (see [KW98]). Three examples of such systems are:

x'[t]-(405/1000) x[t]+(81/100) x[t]  y[t]=0
y'[t] + (15/10) y[t] - (125/1000) x[t]  y[t] = 0
x[0]=6.7   and    y[0]=0.5

Systems such as   (2)  are known to have periodic solutions for a wide range of initial conditions as is easily seen from phase plane analysis. We recommend that this analysis be done before any numerical solutions are generated.  System  (2)  describes a predator-prey model where populations are measured in thousands and time in years. The solution of  (2)  describes (quite accurately according to historical data) an interaction between lions and wildebeest in the Kruger National Park, South Africa (see  [FG99]).  We note that this system cannot be manipulated to represent a second-order ODE.  High-precision approximate closed-form solutions for this model were discussed in  [J01], where the idea to produce this notebook originated.

x'[t] - y[t] = 0
y'[t] + x[t] - ((((x[t]))^3)/6)=(1/3)Cos[0.162615764905768024906 t]
x[0] = 0     and   y[0] = 0

System   (3)  describes a Forced Soft Spring ODE  written as a 2×2 system (see  [Z95]). This system has a periodic solution.

x'[t] - y[t] = 0
y'[t]+ ((((x[t]))^2) - 1) x'[t]+x[t]=0
x[0]=2.008619860874843136422 and   y[0]=0

System  (4)  represents Van der Pol's equation (see  [Z95]). This system has a unique, isolated periodic solution and all solutions converge to the periodic one as t ->∞.  We have chosen the initial conditions so that the trajectory is already close to the limit cycle.

A first approach

Save this notebook under some convenient name so that any part which you inadvertently erase can be readily restored.  Values that you must interactively supply, either by taking an educated guess or by reading off graphs, are indicated by a question mark  "?". 

 Be warned,  DO NOT EVALUATE THE NOTEBOOK, because there are values that are declared as "?"  which you will have to fill in interactively.  If these values are not present, the evaluation will seem to freeze. You will then have to wait for the kernel to process the "garbage" and screen print its error messages before it unfreezes. This might take a while! 
 
 To avoid this situation, evaluate each cell as you interactively progress through the notebook.

For instance, keeping  (1)  in mind,  we will be starting with a routine of the form:

Clear[T, α, β, WP, f, g, res1, res2, sol];
T = ?;
α = ?;
β = ?;
WP = 20;
f[x, y, t] = ?;
g[x, y, t] = ?;
res1[x, y, t]=x'[t]-f[x, y, t];
res2[x, y, t]=y'[t]-g[x, y, t];
sol = NDSolve[{res1[x, y, t] == 0,  
                  res2[x, y, t] == 0,  
                  x[0] == SetAccuracy[α, WP],  
                  y[0] == SetAccuracy[β, WP]}, {x, y}, {t, 0,  
                  4  T}, WorkingPrecision -> WP,  
                MaxSteps -> 5000][[1]]
MaxMemoryUsed[ ]
ParametricPlot[{x[t] /. sol,  
                y[t] /. sol}, {t, 0, T}, PlotPoints -> 100,  
              PlotLabel -> "Phase Plane", PlotPoints -> 100];
Plot[{x[t] /. sol}, {t, 0, T},  
              AxesLabel -> {t, x[t]}, PlotPoints -> 100];
Plot[(y[t] /. sol), {t, 0, T},  
              AxesLabel -> {t, y[t]}, PlotPoints -> 100];

The reader must supply appropriate values and functions. The system the reader is trying to analyse might require different settings for the NDSolve routine such as "Method->Gear"  (for a stiff system), et cetera.

In order to illustrate how the notebook works, we will analyse system   (3).   We leave the analyses of system  (2)  and  (4)  to the reader.  First we will evaluate the NDSolve and graphing routine below which is obtained from  (3)  and  (5).  Because we know that systems  (2),  (3)  and  (4)  have periodic solutions, with period  P*  say,  we expect to be able to estimate  P*  from graphs and to see a closed orbit in the phase plane for  0 ≤  t ≤  P* ([KW98] Method 2).    We have already done an analysis of systems  (2),  (3)  and  (4)  using  a perturbation method (see [KW98] Method 4)  and there is evidence that an initial WorkingPrecision   WP  of  20  is satisfactory for these systems and that the NDSolve routine is working properly for  time  t ≤ 100.   Be warned, if one is analyzing a system other than (2), (3)  or  (4) ,  first do a similar analysis to determine a suitable initial  WP  and time interval.   Further evidence that the  NDSolve algorithm is working properly for systems  (2), (3)  and  (4)   is given in this notebook - a move to higher precision causes the "residuals" to decrease ([KW98] Method 6). Hence we initially set the WorkingPrecision  WP at 20  and the initial maximum "t" value T  is set at some arbitrary positive value which we "tweak" as follows:  If the "Phase Plane"  graph below does not show a closed orbit, choose a larger  T.  Estimate the smallest value for T which is slightly larger  than that which you guess to be the period  P*  of the solution, using the graphs of either x[t] or y[t]  (to get the coordinates, move around in the graphic holding down the Command or Control key (or  equivalent)).  Make sure that a closed orbit is obtained with this vale of  T. 

We worked on a PC with dual 1GHz Pentium III  processors and  2 GB RAM.  Depending on one's computing power one can try   WP>20.

We evaluated the following phase plot and graphs for system   (3)  for   T ≤ 50  in less than a minute:

Determining the approximate period P

Consider that the exact period  P*   is a root of the graphs   G[t] = x[t] - α   and   H[t] = y[t] - β.   Hence evaluate the following commands and examine a graph of  G[t]   and   H[t]   on the same set of axes:

G[t_] := (x[t] /. sol) -  
        SetPrecision[α, WP]
H[
        t_] := (y[t] /. sol) -  
        SetPrecision[β, WP]
Plot[{G[t], H[t]},{t, 0, T},
AxesLabel->{"Time","G[t] is in RED and H[t] is in GREEN"},
PlotPoints -> 100,PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}];

The period  P*  will be the (smallest) value on the  t-axis on the graph above at the point where the curves   G[t],  H[t]  and the t-axis simultaneously intersect and the shape of the curves begins to repeat itself for the first time.  With an eye on the FindRoot routine, we need to determine two numbers   a = ?   and b = ?  close to  P*  satisfying:

(i)    a ≤\(P\^*\)≤b
(ii)  either curve only intercepts the
          t-axis once between these two points

Now enlarge the graph above by clicking and dragging if necessary and read off suitable values of   a   and  b   (to get the coordinates, move around in the graphic holding down the Command or Control key (or  equivalent).  For system  (3)  we have chosen  a  =  38  and  b =39).   Now evaluate:

a = ?
b=?

We now produce two approximate values P1  and  P2  for  P*  (corresponding to the root of  G[t]  and  H[t]  respectively)  by evaluating the following commands:

P1=FindRoot[Evaluate[G[t] == 0],{t, a, a, b}, MaxIterations -> 30,WorkingPrecision -> WP][[1]][[2]]

P2=FindRoot[Evaluate[H[t] == 0], {t, a, a, b},   
              MaxIterations -> 30,  
              WorkingPrecision -> WP][[1]][[2]]

If the FindRoot routine immediately above converges to either of  P1  or  P2   without producing a message and it appears to be what you expect,  skip the next  Plot  evaluation (even if FindRoot  gives messages about the other value).   However, if  both  P1  and  P2  are obviously incorrect or the routine produces the message "Reached the point {c} which is outside the region  {{a, b}}" for both  P1  and  P2,  then zoom in on the point  P  using the next  Plot  command. Choose more appropriate values for a  and  b,  evaluate them as indicated above, re-evaluate the commands for  P1  and  P2  and check to see that the problem has been corrected for at least one of the values  P1  or  P2.

Plot[{G[t], H[t]}, {t, a, b},  
        AxesLabel -> {"Time", "G[t] is in RED and H[t] is in GREEN\
"}, PlotPoints -> 100,  
        PlotStyle -> {RGBColor[1, 0, 0],  
            RGBColor[0, 1, 0]}];

Recall that if  x[t]  and  y[t]  have period  P*,  then α = x[P]   and  β = y[P].  Hence, if there is a choice to be made, which ever one of  P1  or  P2  gives the best approximations  of  α  and  β,  that one should be chosen as the approximate  period  P.  To determine whether to use  P1  or  P2  for  P,  use equation  (7)  below which chooses the larger of the precision indicators  δ1  or  δ2  respectively, as evaluated next:

Clear[δ1, δ2]
α1 =  
        SetPrecision[x[P1] /. sol, 2*WP];
β1 =  
        SetPrecision[y[P1] /. sol, 2*WP];
α2 =  
        SetPrecision[x[P2] /. sol, 2*WP];
β2 =  
        SetPrecision[y[P2] /. sol, 2*WP];
δ1 ->  
      Min[IntegerPart[-Log[10,  
              Abs[SetPrecision[α, 2*WP] - α1]]],  
        IntegerPart[-Log[10,  
              Abs[SetPrecision[β,  
                    2*WP] - β1]]]]
δ1 = %[[2]];
δ2 ->  
      Min[IntegerPart[-Log[10,  
              Abs[SetPrecision[α, 2*WP] - α2]]],  
        IntegerPart[-Log[10,  
              Abs[SetPrecision[β,  
                    2*WP] - β2]]]]
δ2 = %[[2]];

We take note of a precision indicator  δ3  where:

δ3=Max[δ1, δ2]

Evaluate the next input to place   δ3   into the computer's memory:

δ3=Max[δ1, δ2]

If  δ3  is not close to WP-10,  then the NDSolve routine is not working satisfactorily with the settings chosen and it is pointless to proceed with further evaluation of this notebook.

To determine the approximate period  P,  we consider 

 P =?
 
 where  we replace the question mark  ?  with  P1  or  P2  as chosen according to equation  (7)   above.  Hence evaluate:

P = ?

Save the notebook now before proceeding.

Expected precision loss of the solution at  t = P

We now compute a sequence of solutions  solution[n],  for   0 ≤ t ≤ 2T  with  WorkingPrecision ascending in steps of 2 with  n = 0,  2,  4,  ... , k  starting with  solution[0] = sol.  As in KW[96], assuming there is convergence of these solutions towards the exact solution as the WorkingPrecision increases,  we produce graphics which indicate how the AccuracyGoal  varies with time for solution[n],  n = 0,  2,  4,  6.   The "absolute error"  in each of the other solutions as compared to the most accurate solution,  solution[k],  is computed and the "dynamic" accuracy goal,  dynamicAG  is taken to be the negative base 10 logarithm of the absolute error.   Evaluate the following code to obtain two plots, one for  x[t]  and one for  y[t],  that should give an indication of how the AccuracyGoal of each solution  n = 0, 2, 4, 6  varies across the interval  0 ≤ t ≤ P and in particular near  t = P  for  k ≥ 10.  This evaluation took us about 10 minutes with  k = 16.  One might need to lower this to  k = 10  (but not lower than this value!) if the computer runs out of memory.  Ignore any error messages about "machine-sized real numbers".

We now give a graphical analysis similar to that done in [KW96] using a table (Table1).  For each curve, the initial  dynamicAG  should be approximately the  AccuracyGoal  (= WorkingPrecision -10)  of the solution in question.  The curves, cycling through red, green and blue colours, should show evidence (especially those representing solutions with higher AccuracyGoal)  of the convergence of the lower-precision solutions towards the highest-precision solution.  Visualise this evidence as follows:  Drag the plots to enlarge them. Now pick any curve on a plot and imagine the best straight line of fit "by sight" through it (see [FJ01]) and colour it the same colour as the curve. So we have at the bottom an imaginary red line, then, moving up, an imaginary  green  line, an imaginary blue line and finally, at the top, an imaginary red line. For evidence of convergence, each of the lines should have a negative slope, but this slope should become shallower for each line starting from the bottom red line and progressing up towards the top red line.  In particular, the upper curves should show a tendency to have flattened out horizontally.  If convergence is evident,  then the bottom red curve (which starts on the vertical axis in the region of  WP-10)  can be used to estimate, with respect to the (unknown) exact solution, the number of digits of precision lost  for the solution  sol  that we will be using in further calculations.   We would expect the maximum drop in precision over  the interval  0 ≤ t ≤ P   to occur at some local minimum point at time   (Subscript[t, x])  =?   on the  lower red   (Subscript[dynamicAG, x])   graph  and   (Subscript[t, y])  =?   on the lower red  (Subscript[dynamicAG, y])  near or at   t = P   (P  is shown as a point on the time axis). The figures above should confirm this  (to get the coordinates, move around in the graphic holding down the Command or Control key (or  equivalent)).  Hence we determine:

t = (Subscript[t, x])  = ?;
minDynamicx = (Subscript[τ, x])[[1]]

    
t =  (Subscript[t, y])  = ?;
minDynamicy = (Subscript[τ, y])[[1]]
Clear[t, τ]

The expected precision drop for the  solution  sol  over the interval  0 ≤ t ≤ P   is thus given by the following formula and evaluation:

precisionDrop =  
      Max[WP - 10 - IntegerPart[minDynamicx],  
        WP - 10 - IntegerPart[minDynamicy]]

precisionDrop =  
        Max[WP - 10 - IntegerPart[minDynamicx],  
          WP - 10 - IntegerPart[minDynamicy]]

Save the notebook now before proceeding.

Truncated Fourier Series for x[t] and y[t]

Working with the solution sol we now use the mean value theorem to determine Fourier Coefficients for both  x[t]/.sol  and  y[t]/.sol.   The "NIntegrate" routine can be used, but takes longer, produces many error messages and does not produce significantly better results for our examples. If the evaluation of input  (5)  did not produce any messages like "Maximum number of 4000 steps reached at the point   t  == ...", then we choose  S =1500  steps in the mean value theorem.  Usually, using  S = 2000  steps does not produce significantly higher precision as can be seen by experimentation, but does take significantly longer.  We expect these Fourier Series to converge quite rapidly to approximations of  x[t]  and  y[t]  with a precision   p  of about   (δ3 - precisionDrop)  digits (equations  (7)  and  (9) ).  So we are going to determine two positive integers  M  and  R,  as small as possible, which we will use to truncate the respective Fourier Series for  x[t]  and  y[t].   To start with, we take   M = 2×WP  (for no reason but that 2×WP  is a starting value which seems to work in most cases) and then, if necessary, we enlarge  M  so that at least the last  4  coefficients in the tables of Fourier Coefficients determined below are   0.  We start off with  x[t].   This evaluation could take more than 10  minutes.

Clear[M, a, b, S];
M = ?;
S = 1500;
ω = ((2  π)/P);
a[n_]  := (( 2)/S)*
          Sum[Evaluate[x[i*P/S] /. sol]*Cos[n*ω*i*P/S], {i, 0,  
              S - 1}];
b[n_]  :=  
        (2/S)*Sum[
            Evaluate[x[i*P/S] /. sol]*Sin[n*ω*i*P/S], {i, 0,  
              S - 1}];

    
Timing[Table[a[n], {n, 0, M}]];
%[[1]]
(Subscript[Γ, a]) =  
      Chop[%%[[2]], 4*(10^(-(WP - 10)))]

    
Timing[Table[b[n], {n, 1, M}]];
%[[1]]
(Subscript[Γ, b]) =  
      Chop[%%[[2]], 4*(10^(-(WP - 10)))]

Similarly for  y[t]  we evaluate:

Clear[R, A, B, S];
R = ?;
S = 1500;
ω = ((2  π)/P);
A[n_]  :=   
        (2/S)*Sum[
            Evaluate[y[i*P/S] /. sol]*Cos[n*ω*i*P/S], {i, 0,  
              S - 1}];
B[n_]  :=  
        (2/S)*Sum[
            Evaluate[y[i*P/S] /. sol]*Sin[n*ω*i*P/S], {i, 0,  
              S - 1}];

    
Timing[Table[A[n], {n, 0, R}]];
%[[1]]
\(\(\(Γ\_A =  
      Chop[%%[\([2]\)], 4*10\^\(-\((WP - 10)\)\)]\)\(
\)
    \)\)
Timing[Table[B[n], {n, 1, R}]];
%[[1]]
(Subscript[Γ, B]) =  
      Chop[%%[[2]], 4*(10^(-(WP - 10)))]

Now define truncated Fourier Series  X[t]  and  Y[t]   for  x[t]   and   y[t]   respectively as follows:

X[t_]=a[0]/2+(∑_(n = 1)^M)(a[n]*Cos[n*ω*t] + b[n]*Sin[n*ω*t])

Y[t_]=A[0]/2+(∑_(n = 1)^R)(A[n]*Cos[n*ω*t]+B[n]*Sin[n*ω*t])

and evaluate them as follows:

Clear[X, Y]
X[t_] := (Subscript[Γ, a])[[1]]/2 + (∑_(n = 1)^M)((Subscript[Γ, a])[[n + 1]]*
                Cos[n*ω*t] + (Subscript[Γ, b])[[n]]*
                Sin[n*ω*t])
Y[t_] := (Subscript[Γ, A])[[1]]/2 + (∑_(n = 1)^R)((Subscript[Γ, A])[[n + 1]]*
                Cos[n*ω*t] + (Subscript[Γ, B])[[n]]*
                Sin[n*ω*t])

Recall that   x[0]/.sol = α    and  y[0]/.sol = β.   We would expect   X[0]   and   Y[0]    to approximate  α  and  β  to about   δ3   digits of precision (equation  (7)).  The following  evaluation should support this observation:

Testing the approximate solutions

If we graph  the solutions  x[t]/.sol  and  y[t]/.sol   respectively against  X[t]   and  Y[t]  over the interval   0 ≤ t ≤  P we should obtain immediate visual confirmation that we are on the right track. If there is a major difference between the respective functions, two red and green curves will appear. Otherwise, a green curve, laid  perfectly  (to the naked eye)  over the red curve, will appear.

Plot[{(x[t] /. sol), X[t]}, {t, 0, P},  
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]},  
        PlotPoints -> 100];

    
Plot[{(y[t] /. sol), Y[t]}, {t, 0, P},  
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]},  
        PlotPoints -> 100];

In order to gauge how "close"  our truncated  Fourier Series   X[t]  and   Y[t]  are to  x[t]/.sol  and   y[t]/.sol  respectively  for  0 ≤ t ≤ P ,  we define the  "absolute error"    ϵ   in both  X[t]  and   Y[t]  respectively  as follows:

(Subscript[ϵ, 1])[t_] =  
        Abs[(x[t] /. sol) -  
            X[t]];
(Subscript[ϵ, 2])[t_] =  
        Abs[(y[t] /. sol) - Y[t]];

Choosing a suitable m = ?   starting with  m = 0,±1,±2,⋯, we now attempt to plot the functions  |(Subscript[ϵ, 1])|   and   |(Subscript[ϵ, 2])|   below a minimal error bar  (10^(-(Δ - m)))   with

Δ = Min[δ3,δ4]

(equation  (7)  and input  (10)).  The values   α4  and  β4,  given immediately below the following graphs, are the truncated Fourier Series approximations of the initial  conditions  α  and  β   rounded to  Δ- m  digits respectively.  For comparison, the initial  conditions  α  and  β  are given as   α5   and  β5  respectively.

m = ?;

    
Δ =  
        Min[δ3, δ4];

    
N[(10^(-(Δ - m)))]

    
Plot[{(Subscript[ϵ, 1])[t], (10^(-(Δ - m)))}, {t, 0,  
          P}, PlotRange -> {{0, P}, {0,  
              1.25*(10^(-(Δ - m)))}},  
        PlotStyle -> {RGBColor[1, 0, 0],  
            RGBColor[0, 1, 0]}];

    
Plot[{(Subscript[ϵ, 2])[t], (10^(-(Δ - m)))}, {t, 0,  
          P}, PlotRange -> {{0, P}, {0,  
              1.25*(10^(-(Δ - m)))}},  
        PlotStyle -> {RGBColor[1, 0, 0],  
            RGBColor[0, 1, 0]}];

    
Clear[α4, β4, α5, β5]

    
α4 ->  
      SetPrecision[X[0], Δ - m]
α5 -> α

    
β4 ->  
      SetPrecision[Y[0], Δ - m]
β5 -> β

This evidence suggests that these truncated Fourier Series represent closed-form solutions to the system in question over the interval   0≤  t ≤ P  to a precision of    p   digits, where:

p = Δ - m - precisionDrop

p  =  Δ  -  m  -  precisionDrop

Recall that a (vector) function  is a solution of an IVP if it satisfies both the ODE part and the initial conditions. With reference to   (1)  and  (5)   we are going to examine how close the absolute values of the two residual functions 

res1[X, Y, t]   and   res2[X, Y, t] 

are to the zero function, even though, as opposed to the exact solution,  this is generally not an indication of the validity of an approximate solution (see [KW98] Method 6) . 

Choosing a suitable  q  =  ?,  starting with  q = 0,±1,±2,⋯, we now attempt to plot the functions  |res1|   and   |res2|   below a minimal error bar  (10^(-(Δ - q))) :

q = ?;
N[(10^(-(Δ - q)))]

Plot[{Abs[res1[X, Y, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                  1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  X[t]",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, PlotPoints -> 100];

Plot[{Abs[res2[X, Y, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  Y[t]",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0],  
              RGBColor[0, 1, 0]}, PlotPoints -> 100];

It is interesting to observe that the absolute values of the residuals of the NDSolve solutions  x[t]/.sol  and  y[t]/.sol  are  usually larger than  (10^(-(Δ - q))).   In order to see this, evaluate the next input.

N[(10^(-(Δ - q)))]

Plot[{Abs[res1[x /. sol, y /. sol, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                  1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  x[t]/.sol",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, PlotPoints -> 100];

Plot[{Abs[res2[x /. sol, y /. sol, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  y[t]/.sol",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0],  
              RGBColor[0, 1, 0]}, PlotPoints -> 100];

We also note that a move to higher precision causes the NDSolve residuals to decrease, further evidence that the NDSolve routine is working properly for the system being analysed ([KW98] Method 6) as can be seen by evaluating the next plot of the residuals for  solution[4]   (see input (8)):

solution[4];
sol4 = %;

N[(10^(-(Δ - q)))]

Plot[{Abs[res1[x /. sol4, y /. sol4, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                  1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  x[t]/.solution[4]",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, PlotPoints -> 100];

Plot[{Abs[res2[x /. sol4, y /. sol4, t]],(10^(-(Δ - q)))},{t, 0, P},PlotRange -> {{0, P}, {0,  
                1.25*(10^(-(Δ - q)))}},PlotLabel -> "Residual for  y[t]/solution[4]",PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0],  
              RGBColor[0, 1, 0]}, PlotPoints -> 100];

Conclusion

Keeping equation  (9)  in mind, assuming that the NDSolve solution  sol  approximates the exact solution up to  (WP-10-precisionDrop)   digits over the interval  0 ≤  t  ≤  P,   we have produced approximate solutions in the form of truncated Fourier Series X[t]  and  Y[t]  with a precision of  p   digits  (equation  (12))  over the interval  0 ≤  t  ≤  P.  If,  for instance, X[t]  and  Y[t]  represent populations measured in thousands, then  p ≥ 4  should be sufficient  for predictions over the interval   0 ≤  t  ≤  P .   If  p  is too low for the physical situation being modelled, then the notebook should be re-evaluated, starting with an appropriately larger WorkingPrecision  WP.   Because the truncated Fourier Series  X[t]  and  Y[t]  periodically extend the solution  sol  for   t  ≥  P  and because we expect  Mathematica (or a hand-held scientific calculator that can handle  p  digits  of precision) to maintain precision for a truncated Fourier Series,  we could also expect  X[t]  and  Y[t]   to approximate the exact solution to   p  digits for t  ≥  4P,   for a while.  Indeed, if  P*  is the period of the exact solution,  unless the approximate period  P = P*,  one cannot expect the truncated Fourier Series to approximate the true solution forever. The following observation should make the last remark clearer.  If   |P*-P| ≈ (10^(-r))   say,  then after  10P*  time periods the graph of  X[t]  will have translated horizontally to the left or the right through a distance of approximately (10^(-r + 1))  from  the exact solution  x[t],  and similarly for  Y[t].  

One final observation:   By re-examining the graphs of the absolute errors  (Subscript[ϵ, 1])  and  (Subscript[ϵ, 2])  plotted out to  t = 4T,  one may observe how the NDSolve solutions  x[t]/.sol  and  y[t]/.sol  diverge from  X[t]  and  Y[t]  respectively over the interval  0 ≤ t ≤  4T  (the period  P  is shown as a point on the time axis).  Indeed, evaluate the following code:

Plot[Abs[(Subscript[ϵ, 1])[t]], {t, 0, 4  T},  
        PlotStyle -> {RGBColor[1, 0, 0]}, PlotPoints -> 100,  
        Epilog -> {AbsolutePointSize[10],  
            Point[{P, 0}]}];

    
Plot[Abs[(Subscript[ϵ, 2])[t]], {t, 0, 4 T},  
        PlotStyle -> {RGBColor[1, 0, 0]}, PlotPoints -> 100,  
        Epilog -> {AbsolutePointSize[10], Point[{P, 0}]}];

All calculations and visualizations have been carried out in Mathematica 4.1.

Bibliography

[FJ01]:  Temple H. Fay and Stephan V. Joubert, Least squares and energy, Mathematics and Computer Education, 35(2), (2001), 114-115. 
 [J01]: Stephan V. Joubert, High Precision Approximate Closed Form Solutions For a Predator-Prey Model, Proceedings of the WarthogDelta Conference, Kruger National Park, South Africa, 2001.
[FG99] : Temple H. Fay and Johanna C. Greeff, Lions and wildebeest: A predator-prey model, Mathematics and Computer Education, 33(2) (1999), 106-119.
 [KW98]:  Rob Knapp and Stan Wagon, MATHEMATICA IN ACTION: Check Your Answers . . . But How? ,  Mathematica in Education and Research, 7(4), (1998), 76-85.
 [KW96]: Rob Knapp and Stan Wagon, Orbits worth betting on, C.ODE.E, Winter issue, (1996), 8-13.
 [Z95]: Dennis G. Zill, Differential Equations with Computer Lab Experiments, PWS Publishing Company, Boston, 1995.