# Stiff3 - adaptive solver for stiff systems of differential equations

stiff3 is a semi-implicit Runge Kutta solver with adaptive time-stepping originally from the book by Villadsen & Michelsen, Solutions of Differential Equation Models by Polynomial Approximation, Prentice-Hall, 1978.

The solver is now available as an fpm package and uses LAPACK for the matrix factorization and back-substitution in the semi-implicit time-stepping.

The solver is designed for autonomous problems of the form

\frac{dy}{dx} = f(y)

In case of a non-autonomous problem where x appears explicitly, i.e. y' = f(x,y), it suffices to switch to a new integration variable t and introduce the additional ODE

\frac{dx}{dt} = 1

The repository provides an example for the Van der Pol oscillator. Here are some plots of the solutions and phase cycle.

There is one remaining issue with the stiff3 solver where I’m seeking help. It has to do with a section of code that sets to zero some elements of the Jacobian matrix. I believe this section is just an outdated practice that has to do with legacy usage patterns. I’ve posted it as a question on SciComp StackExchange, but so far it is still unanswered.

6 Likes

Nice! Thanks for creating the library. We need more of such libraries.

You could code the do-loops in question as:

where ( abs(df) > df_tol )
df = -h * a * df
elsewhere
df = 0.0_wp
endwhere

forall (i=1:n)
df(i,i) = df(i,i) + 1.0_wp
end forall


(If I remember the forall syntax properly)

Whether you actually want to do it that way, depends on your appreciation of these constructs

Indeed, that is more elegant. However, I think the zero values should be set explicitly in the user-provided Jacobian routine, and there is no need for this chunk of code at all. It doesn’t make sense to me to hardcode a df_tol cutoff in the Jacobian matrix.

Is LAPACK supported in FPM now ?

Not exactly. But you can add it as a build dependency:

[build]

As long as you have a system-installed LAPACK library that the linker is able to locate with flag -llapack, then things should just work.