From: Panagiotis Papasotiriou <p.j.papasot@gm...>  20140420 13:53:44

The best way to deal with stiff ODEs would be a multistep method able to switch between nonstiff and stiff ODE solvers automatically. As already mentioned by others here, ODEPACK <http://www.netlib.org/odepack>; is probably the way to go. It is actually the best package for that purpose, to the best of my knowledge. ODEPACK's solver named LSODAR (and its double precision variant, DLSODAR) is not only able to switch between "stiff" and "nonstiff" mode automatically, but also has builtin root finding capabilities: the user can optionally supply a set of constraint functions (involving both dependent and independent variables), and (D)LSODAR will stop integration whenever a root of such a constraint function is found. This feature is very useful when the integration domain is not known in advance. As a Fortran programmer, I used DLSODAR a lot, and I have yet to see it failing solving a initialvalue problem. Packages like Scilab or the R language use ODEPACK internally to solve initialvalue problems, and I have reasons to believe commercial packages like Matlab or NAG libraries do the same (the few functions I had access to, used in such commercial products, were actually open source routines renamed, verified by comparing them with code found in netlib.org.) Based on the above, I believe that being able to do use DLSODAR in Maxima would be the best way to deal with both stiff and nonstiff ODEs within Maxima. I don't expect it to be as fast as using a native executable, but I do expect the same robustness. DLSODAR, together with its (stripped) ODEPACK and SLATEC dependencies, is more than 4500 lines of code, but in my opinion it is well worth the effort. Besides, I' have read here that ODEPACK is already translated somehow to Lispreadable code. Although both my editor and CAS of choice (Emacs and Maxima) are written in Lisp, I am not familiar with that language, sadly. The few times I gave it a try, I ended up thinking "why should I do it that way when I could use Fortran 95/2003 instead?" (plus Emacs' Lisp mode, "Slime", never worked well in my system.) Thus I am not able to provide any help in the Lisp domain, but I am able to provide a Maxima interface function to DLSODAR. ODEPACK itself is written in old Fortran 77 spaghetti code (full of GOTOs,) but I've already written a Fortran 95 interface<https://sites.google.com/site/pjpapasot/fortran/libraries/dlsodar_f95>to DLSODAR, which provides all the original functionality but automatically sets most of DLSODAR's many arguments, and in general makes things easier. In general, translating that interface to Maxima for easy use doesn't seem to be that hard, provided a Maximareadable version of ODEPACK does exist. 20140417 23:39 GMT+03:00 Edwin Woollett <woollett@...>: > On April 16, 2014, Helmut Jarausch wrote: > > > On 04/09/2014 10:28:11 PM, Edwin Woollett wrote: > >> > >> feature request: stiff o.d.e. solvers for initial value problems > >> in Maxima. > > > To teach myself coding in Maxima, I have translated a code which I have > > written in Scilab/Matlab > > for teaching purposes several years, ago. > > > It's an extrapolated, linearlyimplicit Euler method. I have written a > > simplified version > > of Ernst Hairer's SEulex FORTRAN code. > > > It features variable step size, variable order and dense output. > > It can switch from explicit to implict and vice versa, but not > > automatically. > > It can solve DAEs of index 1, as well. > > > The code seems to work, though a bit slow. > > > I appreciate any comments, > > > See the examples and the core code at > > > http://www.igpm.rwthaachen.de/jarausch/Maxima > >  > > This code needs a userfriendly interface which has default values > for the many options. > > The user should be able to get back a solution list (as in rk or rkf45) > after providing the same kinds of input rkf45 needs to start work. > > In other words, for 1 o.d.e. > > eulix(dydt, y, y0, [t, t0, tf]) > > and for 2 o.d.e.s, > > eulix([dy1dt, dy2dt], [y1,y2],[y10,y20], [t, t0, tf]) > > So your "mass matrix" would automatically be a unit matrix, and the user > should not be concerned with matrices. > > As far as setting defaults for optional args, > see the code for rkf45 which uses assoc. > for example: >  > /* Set optional arguments */ > atol:assoc('absolute_tolerance,options,1e6), > save_steps:assoc('full_solution,options,true), > maxit:assoc('max_iterations,options,10000), > show_report:assoc('report,options,false), > > etc, etc. >  > Express your solved examples in ordinary language, such as > solve dy/dt = t*y, with y(0) = 2, over range [t,0,5], so the > file reader does not have to trace through your code to > understand how to use it. > > Start with the simplest goal of getting the code to work > without bells and whistles at first. Try to let Maxima > do any hard matrix work with core matrix methods > rather than translating fortran matrix routines. > > Consider breaking up the program into smaller chunks which > can be separately debugged. > > Thanks for your efforts on this project. > > Ted Woollett > > > > > > >  > Learn Graph Databases  Download FREE O'Reilly Book > "Graph Databases" is the definitive new guide to graph databases and their > applications. Written by three acclaimed leaders in the field, > this first edition is now available. Download your free book today! > http://p.sf.net/sfu/NeoTech > _______________________________________________ > Maximadiscuss mailing list > Maximadiscuss@... > https://lists.sourceforge.net/lists/listinfo/maximadiscuss > 