# Is it worth switching to timesteppers provided by PETSc if I can't write down a Jacobian for my problem? Case study with “the amoeba” toy problem

I am considering using petsc4py instead of scipy.integrate.odeint (which is a wrapper for Fortran solvers) for a problem involving the solution of a system of ODEs. The problem has the potential to be stiff. Writing down its Jacobian is very hard.

So far, I have been able to produce reasonable speed gains by writing the RHS functions in "something like C" (using either numba or Cython). I'd like to get even more performance out, hence my consideration of PETSc.

# Introducing a toy version of the problem: The Amoeba

Consider the following figure, which defines the basics of an amoeba:

So, an amoeba consists of $N$ nodes and edges, constrained to a 2D space. The position of the nodes describing the amoeba evolve with time.

Let the position of the $n$th ($n \in \{0, 1, ..., N-1\}$) node be described a position vector, $\textbf{x}_n$:

$$\textbf{x}_n = \left[\begin{array}{cc} x_{0,n} \\ x_{1,n} \end{array}\right]$$

For the purposes of simplification, assume that the amoeba lives in a world with first order equations of motion, so given an applied force $\textbf{F}$ on a node, the nodal velocity $\dot{\textbf{x}}$ is produced:

$$\textbf{F} = c \dot{\textbf{x}}$$

where $c$ is a constant I pick.

We are very interested in tracking the position of each node,so immediately we produce $2n$ ODEs (one ODE for each spatial component) of interest (assume that edges are elastic springs):

$$\frac{\textrm{d}\textbf{p}_n}{dt} = \text{elastic edge forces} + \textbf{F}_{n, active}$$

$\textbf{F}_{n, active}$ is the force generated by the amoeba at node $n$. One can imagine a wide range of rules that could be devised to determine how $\textbf{F}_{n, active}$ is determined, and I suggest that we consider the following rules:

1) There are red balls and blue balls inside the amoeba. ^ 2) Both red and blue balls have "switched-on" (symbolized $r^*$/$b^*$) and "switched-off" (symbolized $r$/$b$) state.

3) Both red and blue balls can be found either in the body of the amoeba where they are always switched off, or at the nodes of the amoeba, where they can be switched on or switched off.

4) The total number of red balls regardless of location of status is a constant $R$ and similarly for blue balls, a constant $B$

5) When red or blue balls are in the body, they are free to jump onto any node they want, but when they are on nodes, some of them might be induced to move to neighbouring nodes due to some sort of "diffusion"

6) Switched on red balls at a node cause a force that attempts to reduce the area of the amoeba

7) Switched on blue balls at a node cause a force that attempts to increase the area of the amoeba

8) Switched on red balls tend to switch off blue balls that are switched on, and vice versa, switched on blue balls tend to switch off red balls that are switched on

9) ...

# Okay, you get my point. I can come up with a model that is fairly convoluted in terms of how the various variables of interest interact.

Due to the large number of equations involved, it is already tedious to think about writing downa Jacobian. However, since we have access to a computer, it may be that we can write functions governing a particular interaction that do not have neat analytical forms (let alone whether or not their derivatives have neat analytical forms), so we might have a mess of piecewise functions needed to approximate them if we were to go about still trying to produce a Jacobian...

# Things get even worse if my problem involves not just one amoeba, but multiple amoebae (which is when things ARE interesting), with new interaction rules enforcing that amoebae cannot violate each others' areas, and so on.

All the toy examples I see of PETSc time stepping problems have Jacobians defined, so I wonder if I would even get a speed gain going from switching to it, if perhaps one of the reasons why I have a high computational cost is due to not being able to provide a Jacobian function?

• Have you considered to calculate your Jacobian by automatically differentiating the (C code for the) right-hand-side of your ODE system? – GoHokies Dec 9 '15 at 13:54
• As a followup to GoHokies's comment, I've used DAE Tools to set up systems of ODE's (or DAE's) in python and let it take care of automatic differentiation (and some other details) for me. – muon Jul 29 '16 at 13:57