BDF2a: An adaptive-stepsize Stiff ODE solver in Maxima

Since 2011, Maxima has included the user-contributed numerical ODE solver rkf45() created by Panagiotis Papasotiriou.  This implementation of the fourth (and fifth) order Runge-Kutta-Fehlberg embedded method features adaptive timestep selection and a nicely optimized function evaluation to make it run pretty fast in Maxima.  As an explicit RK method, it is suitable for non-stiff equations.  To use it in Maxima, you must first load(rkf45).

For years in my differential equations class, we’ve test-driven the various numerical ODE methods in MATLAB and explored their behavior on stiff systems.  This year we’ve switched from MATLAB to Maxima, and I’ve written a stiff solver for that purpose.  It’s an adaptive stepsize  implementation of the second-order Backward Difference Formula.  The BDF.mac file can be found here.  The package contains a fixed-stepsize implementation BDF2()  and the adaptive stepsize BDF2a().

I’ve tried to make the calling sequence and options as close to rkf45() as possible.  I have not made any attempt to optimize the performance of the method.  BDF methods are implicit, meaning that a nonlinear equation (or system of equations) must be solved at each step.  For that I use the user-contributed Maxima function mnewton().

Additionally, note that I included my homemade help utility so that once BDF.mac is loaded, help(BDF2) and help(BDF2a) display help lines in the console window. For simplicity in comparing the methods, the load(rkf45) command is included in BDF.mac, as is an rkf45 help file for use with help().   Finally, I included in the package my little solution plotting utility wxtimeplot() which takes the “list of lists” output of BDF2(), BDF2a() and rkf45() and plots the dependent variables vs the independent variable along with a log-scale plot of the stepsize.

Below is an example that illustrates how the the BDF method performs relative to rkf45 for an increasingly stiff problem.

The equation is

 \dot{y} = -\lambda (y-\sin (10t)-t)+10 \cos (10t)+1 .    

For small values of  \lambda , the explicit rkf method’s adaptive stepsizes result in much more efficient solutions for a given accuracy level.  However for large \lambda the implicit BDF method is comparatively much more efficient, in terms of number of steps, and  number of function evaluations, although in its current implementation, the BDF method takes more time.

Here’s the results for \lambda=50:

stiff50

stiff50plot

But for \lambda=5000, the performance situation is reversed, with the BDF formula essentially immune from the effects of extra stiffness, while the rkf method requires an order of magnitude more steps.

stiff5000

stiff5000plot

References:

E. Alberdi Celaya , J. J. Anza Aguirrezabala , and P. Chatzipantelidis.  “Implementation of an adaptive BDF2 formula and comparison with the MATLAB Ode15s”,  Procedia Computer Science,  December 2014

 

 

 

Advertisements

2 thoughts on “BDF2a: An adaptive-stepsize Stiff ODE solver in Maxima”

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s