Finite difference method converging to wrong steady state (python)

Question

I'm trying to recreate numerical results from an article. In particular, I want to numerically solve the following PDE system:

In this system, Fn represents a prey population density, Cn a predator population density, and En the energy of the predators. The (3.57) equations model how these different populations interact when time t lies between two consecutive natural numbers, i.e., n<t<n+1 for some natural number n. The (3.4) equations govern how the populations change when time is a natural number, i.e., when t=n for some natural number n.

When the authors numerically solve the non-spatial system, i.e., with all the diffusion terms in the (3.57) equations dropped, they produce the following results (the parameters they use are at the bottom of the image):

When I numerically solve the system in Python, I get the following results:

My results look very similar, however, the peaks of my solutions are always around twice as high as the solutions from the article. This discrepancy doesn't disappear by taking smaller time steps so I'm thinking something is wrong with how I'm implementing the forward Euler method in Python. Any help would be greatly appreciated.

import matplotlib.pyplot as plt

ρ=0.55
r1=10
θ1=0.17
θ2=0.1
m1=2.1
β=7
ω=0.01
δ1=0.05

h=0.01 #step size
T=15000 #number of steps

D=[] #domain

for i in range(0,T+1):
    D.append(i*h)

F=0.05 #prey initial condition
C=2.5 #predator initial condition
E=0 #energy initial condition

Fval=[F] #list that will store all calculated prey values
Cval=[C] #list that will store all calculated predator values
Eval=[E] #list that will store all calculated energy values

for i in range(1,T+1):
    if (i*h)%1==0: #between year dynamics
        Fnew=ρ*F
        Cnew=E+C
        Enew=0
    else: #within year dynamics
        Fnew=(1+h*r1*(1-F))*F-((h*F*C)/(ω+θ1*F+θ2*C))
        Cnew=(1-h*m1)*C
        Enew=((h*β*F*C)/(ω+θ1*F+θ2*C))+(1-h*δ1)*E
        
    Fval.append(Fnew)
    Cval.append(Cnew)
    Eval.append(Enew)

    F=Fnew
    C=Cnew
    E=Enew

#plotting
plt.plot(D, Fval,color='red',linestyle='dotted')
plt.plot(D, Cval,color='blue',linestyle='dashdot')
plt.plot(D, Eval,color='black')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$, $C_n(t)$, $E_n(t)$')
plt.legend(['$F_n(t)$', '$C_n(t)$','$E_n(t)$']) 
plt.xlim(0,50)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial homogeneous impulsive solutions')
plt.show()

plt.plot(D, Fval,color='blue')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Cval,color='blue')
plt.xlabel('time t')
plt.ylabel('$C_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Eval,color='blue')
plt.xlabel('time t')
plt.ylabel('$E_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

Answer 1

Well, I tried coding it up independently. Sorry, it's in C++, with graphs from gnuplot. It's also using 4th-order Runge-Kutta.

Basically, ... I agree with YOUR results!

The final state doesn't seem sensitive to initial conditions, so I think the authors were perhaps using different values of the coefficients - perhaps you could try changing them.

Since at least this gives you independent confirmation of your result it might be worth contacting the authors and asking them to check their coefficients.

#include <iostream>
#include <valarray>
using namespace std;

const double rho = 0.55, r1 = 10.0, theta1 = 0.17, theta2 = 0.1, m1 = 2.1, b = 7, omega = 0.01, delta1 = 0.05;

//======================================================================

valarray<double> rk4( double x, const valarray<double> &Y, valarray<double> f( double, const valarray<double> & ), double dx )
{
   valarray<double> dY1, dY2, dY3, dY4;
   dY1 = dx * f( x           , Y             );
   dY2 = dx * f( x + 0.5 * dx, Y + 0.5 * dY1 );
   dY3 = dx * f( x + 0.5 * dx, Y + 0.5 * dY2 );
   dY4 = dx * f( x +       dx, Y +       dY3 );
   return Y + ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
}

//======================================================================

valarray<double> deriv( double t, const valarray<double> &Y )
{
   double F = Y[0], C = Y[1], E = Y[2];
   double dFdt = r1 * ( 1.0 - F ) * F - F * C / ( omega + theta1 * F + theta2 * C );
   double dCdt = -m1 * C;
   double dEdt = b * F * C / ( omega + theta1 * F + theta2 * C ) - delta1 * E;
   return valarray<double>{ dFdt, dCdt, dEdt };
}

//======================================================================

int main()
{
   valarray<double> Y = { 0.05, 2.5, 0.0 };
   const int N1 = 100;   // timesteps per year
   const double dt = 1.0 / N1;
   const int Years = 150.0;
   double t = 0;
   for ( int y = 0; y < Years; y++ )
   {
      // Within a year
      for ( int n = 0; n < N1; n++ )
      {
         cout << t << '\t' << Y[0] << '\t' << Y[1] << '\t' << Y[2] << '\n';
         Y = rk4( t, Y, deriv, dt );
         t += dt;
      }
      // End-of-year reset
      Y[0] = rho * Y[0];
      Y[1] = Y[2] + Y[1];
      Y[2] = 0.0;
   }
}

Graphs

Answer 2

I believe that the nonlinearity is the problem. Forward Euler in the form that you have works great for linear ODEs, but your first ODE is quite heavily nonlinear. (I will call these ODEs, because when you remove the spatial diffusion, these are just couple equations in t only.)

There appears to be nonlinear damping via the F_n^2 term, plus resonant damping in via the F_n C_n term. Both of these will be quite sensitive to change, and too sensitive to smooth out by shrinking the grid size.

My recommendation is to see if the authors do indeed use forward Euler (it's not always suitable!) From experience, I would typically try solving this using RK4 (4th order Runge-Kutta) first.

RK4 is generally quite good at solving nonlinear equations similar to this, provided you use a fine discretization in time.

Good luck!