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.

image

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. :slight_smile:

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 :slight_smile:

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]
link = ["lapack", "blas"]

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. :rocket:

Great stuff @ivanpribec! Nice to see the fpm ecosystem growing with more useful packages!

Regarding your final point, I’ve previously come across thresholding of small matrix elements for iterative methods - I believe it’s to avoid having very small eigenvalues which degrade convergence of the overall method or place severe limits on the region of stability.

As pointed out on the scicomp post, these small values creep in unavoidably due to floating-point arithmetic, but they are essentially zero because they are smaller than the effective numerical precision which gets smaller as matrix condition number increases. (I think I read somewhere that every order of magnitude in the condition number takes off one significant figure in numerical precision.) The value of the cut-off threshold was perhaps set based on some expected condition number.

2 Likes

Very nice, thank you. I suggest submitting stiff3 to GitHub - fortran-lang/fpm-registry: Centralized registry of fpm packages.