It is assumed that the ODEs are given explicitly, so that the system can be written in the form dy/dt = f(t,y), where y is the vector of dependent variables, and t is the independent variable. The Krylov methods are based on Arnoldi's method, the Generalized Minimum Residual method, and the Preconditioned Conjugate Gradient method, with scaling added.
This program solves a stiff ODE system that arises from a system of partial differential equations. The PDE system is a food web population model, with predator-prey interaction and diffusion on the unit square in two dimensions. The dependent variable vector is
1 2 ns
c = (c , c , ..., c )
and the PDEs are as follows:
i i i
dc /dt = d(i)*(c + c ) + f (x,y,c) (i=1,...,ns)
xx yy i
i ns j
f (x,y,c) = c *(b(i) + sum a(i,j)*c )
The number of species is ns = 2*np, with the first np being prey and the last np being predators. The coefficients a(i,j), b(i), d(i) are:
a(i,i) = -a (all i)
a(i,j) = -g (i .le. np, j .gt. np)
a(i,j) = e (i .gt. np, j .le. np)
b(i) = b*(1 + alpha*x*y) (i .le. np)
b(i) = -b*(1 + alpha*x*y) (i .gt. np)
d(i) = dprey (i .le. np)
d(i) = dpred (i .gt. np)
The various scalar parameters are set in subroutine setpar.
The boundary conditions are: normal derivative = 0.
A polynomial in x and y is used to set the initial conditions.
The PDEs are discretized by central differencing on a mx by my mesh.