Lotka–Volterra equations

[beantsja]

How would I write an equation for Lotka-Volterra equation? It is because I am planning to expand this to writing a four coupled PDE that has a source term similar in this form. For example, I am solving for x,y,z,g, and the source terms are like xy, yz and etc...

eqn0 = TransientTerm(var=x) == a*ImplicitSourceTerm(coeff=1.0, var=x)-b*ImplicitSourceTerm(coeff=1.0, var=x)*ImplicitSourceTerm(coeff=1.0, var=y)
eqn1 = TransientTerm(var=y) == d*ImplicitSourceTerm(coeff=1.0, var=x)*ImplicitSourceTerm(coeff=1.0, var=y)-c*ImplicitSourceTerm(coeff=1.0, var=y)

eqn = eqn0 & eqn1

This is what I have, but it is giving me an error called "TermMultiplyError: Must multiply terms by int or float."

[guyer]

You cannot multiply a Term by another Term.

FiPy uses the Finite Volume method to discretize PDEs into linear algebra equations. Intrinsic to linear algebra is that Terms can only be implicit in one variable, so you must choose.

I would write

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * x, var=y)
eqn1 = TransientTerm(var=y) == ImplicitSourceTerm(coeff=d * y, var=x) - ImplicitSourceTerm(coeff=c, var=y)

eqn = eqn0 & eqn1

(The alternative,

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * y, var=x)
eqn1 = TransientTerm(var=y) == ImplicitSourceTerm(coeff=d * x, var=y) - ImplicitSourceTerm(coeff=c, var=y)

eqn = eqn0 & eqn1

would result in two equations that aren't coupled at all. While they would probably solve, they would probably converge more poorly.)

[guyer]

Beyond that, what you have so far are ODEs rather than PDEs. I'd think that FiPy could solve them, but tools like scipy.integrate.solve_ivp() or SUNDIALS would probably do better.

[guyer]

There are a few things going on here:

• Your equations are nonlinear, so they require sweeping to reach convergence. I find that three sweeps is sufficient to reach convergence for your equations.

  48c48,50
  <     eqn.solve(dt=1)
  ---
  >     for sweep in range(3):
  >         res = eqn.sweep(dt=1)

• By setting the initial value to 20, you're storing integers in y. There are reasons to do that, but integers don't work well for solution variables. I'm not sure why you're not getting a warning about that; perhaps because you've coupled with x containing floats.

  9c9
  < y = CellVariable(mesh=m, hasOld=True, value=20)
  ---
  > y = CellVariable(mesh=m, hasOld=True, value=20.)

• With smaller timesteps, the solutions are only (apparently) linear over the same hundred steps. If you, e.g., increase to 1000 steps of 0.1 (for the same total time), then you get fluctuating solutions. With dt=1, I find that the amplitude of y damps over time.
  
  Decreasing the timestep to dt=0.1
  
  or dt=0.01
  
  leads to solutions with more consistent amplitude.
  Smaller timesteps do improve accuracy, but they are quite slow. This is both because FiPy isn't designed for solving ODEs and because you're solving the ODE 100 times (at every grid point).

[guyer]

Your planned equations

$$ \begin{aligned} \frac{\partial a}{\partial t} &= D_a \frac{\partial^2 a}{\partial x^2} + k b a \\ \frac{\partial b}{\partial t} &= D_b \frac{\partial^2 b}{\partial x^2} - k b a \end{aligned} $$

are definitely more suitable for FiPy. These would be written as

eqn0 = TransientTerm(var=a) == DiffusionTerm(coeff=Da, var=a) + ImplicitSourceTerm(coeff=k * a, var=b)
eqn1 = TransientTerm(var=b) == DiffusionTerm(coeff=Db, var=b) - ImplicitSourceTerm(coeff=k * b, var=a)

[beantsja]

Thank you very much.

Lastly, could you explain how you derive the equations?

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * x, var=y)
eqn1 = TransientTerm(var=y) == ImplicitSourceTerm(coeff=d * y, var=x) - ImplicitSourceTerm(coeff=c, var=y)

eqn = eqn0 & eqn1

You put x as a coefficient in one of the ImplicitSourceTerm and y for another. Does it matter which variable I choose to put as a coefficient? In four coupled PDE case, can I put more than one variables in the coefficient? For example,

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * x*z, var=y)
considering it is solving for x, y, z.

[guyer]

You can write it either way, but it matters. ImplicitSourceTerm(coeff=f, var=g) represents $f g$. f can be anything, but the expression can only be implicit (linear) in one variable, g.

By choosing

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * x, var=y)
eqn1 = TransientTerm(var=y) == ImplicitSourceTerm(coeff=d * y, var=x) - ImplicitSourceTerm(coeff=c, var=y)

I chose to represent

$$ \begin{align} x' &= a x - (b x) y \\ y' &= (d y) x - c y \end{align} $$

or

$$ \begin{align} (\partial_t - a) x - (b x) y &= 0 \\ -(d y) x + (\partial_t - c) y &= 0 \end{align} $$

which creates a (very simplified) block matrix

$$ \left[\begin{array}{c|c} \partial_t - a & -b x \\ \hline -d y & \partial_t - c \end{array}\right] \left[\begin{array}{c} x \\ y \end{array}\right] = 0 $$

which represents pretty good coupling.

If you write

eqn0 = TransientTerm(var=x) == ImplicitSourceTerm(coeff=a, var=x) - ImplicitSourceTerm(coeff=b * y, var=x)
eqn1 = TransientTerm(var=y) == ImplicitSourceTerm(coeff=d * x, var=y) - ImplicitSourceTerm(coeff=c, var=y)

then you are saying

$$ \begin{align} x' &= a x - (b y) x \\ y' &= (d x) y - c y \end{align} $$

or

$$ \begin{align} (\partial_t - a - b y) x &= 0 \\ (\partial_t - d x - c) y &= 0 \end{align} $$

which creates a block matrix

$$ \left[\begin{array}{c|c} \partial_t - a - b y & 0 \\ \hline 0 & \partial_t - d x - c \end{array}\right] \left[\begin{array}{c} x \\ y \end{array}\right] = 0 $$

The first equation solves only for x and the second solves only for y. You can "couple" them with &, but they're not really coupled and they won't converge any better than if you swept them sequentially. It'll still work, but it will probably converge more slowly.

Actually, it looks like the convergence is about the same either way, with the residual reaching machine precision in around 20 sweeps. I'd originally identified 3 sweeps convergence when y was still an integer. To be honest, the solution doesn't look any "better" after 20 sweeps. Accuracy (assuming constant amplitude oscillations is "accurate") is only achieved by reducing the time step.

[guyer]

Thinking about this a little more, the fact that the off-diagonal blocks ($-bx$ and $-dy$) are negative probably means that they get put on the right-hand-side vector instead of on the matrix (for stability reasons). This is probably why the two different approaches have similar convergence. Still, I'd do the first one; the second is definitely not coupled, whereas FiPy will try to couple the first if it can.