## Re: [Maxima-discuss] Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems]

 Re: [Maxima-discuss] Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems] From: Edwin Woollett - 2014-04-17 20:39:40 ```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, linearly-implicit 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.rwth-aachen.de/jarausch/Maxima --------------------------------------------------------------------------------- This code needs a user-friendly 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,1e-6), 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 ```

 [Maxima-discuss] feature request: stiff ode solvers for initial value problems From: Edwin Woollett - 2014-04-09 20:28:02 ```feature request: stiff o.d.e. solvers for initial value problems in Maxima. Many advanced numerical algorithms that solve differential equations are available as (open- source) computer codes, written in programming languages like FORTRAN or C, and available in repositories like GAMS (http://gams.nist.gov) or NETLIB (www.netlib.org). An example of what can be done to make this code available for work in modern interactive numerical environments is the work of Karline Soetaert and L.R. Petzold (and others) for the open source R numerical platform. Present ode solvers in the R package deSolve use adaptive step size control, some solvers control the order of the formula adaptively, or switch between different types of methods, depending on the local properties of the ode's to be solved. The R package deSolve includes methods to solve stiff and non-stiff problems, that deal with full, banded, or arbitrarily sparse Jacobians, etc. The implementation includes stiff and nonstiff integration routines based on the ODEPACK FORTRAN codes (Hindmarsh 1983). It also includes fixed and adaptive timestep explicit Runge-Kutta solvers, the Euler method, and the implicit Runge-Kutta method RADAU (Hairer and Wanner 2010). The default integration method, based on the FORTRAN code LSODA is one that switches automatically between stiff and non-stiff systems (Petzold 1983). This is a very robust method, but not necessarily the most efficient solver for one particular problem. See (Soetaert et al. 2010b) for more information about when to use which solver in deSolve. For most cases, the default solver, ode, and the default settings will do. One can optionally trigger specific available methods to override the default used by the deSolve wrapper ode. The explicit Runge-Kutta methods are de novo implementations in C, based on the Butcher tables (Butcher 1987). They comprise simple Runge-Kutta formulae (Euler's method euler, Heun's method rk2, the classical 4th order Runge-Kutta, rk4) and several Runge-Kutta pairs of order 3(2) to order 8(7). The embedded, explicit methods are according to Fehlberg (1967) (rk..f, ode45), Dormand and Prince (1980, 1981) (rk..dp.), Bogacki and Shampine (1989) (rk23bs, ode23) and Cash and Karp (1990) (rk45ck), where ode23 and ode45 are aliases for the popular methods rk23bs resp. rk45dp7. With the following statement all implemented rk methods are shown: > rkMethod() [1] "euler" "rk2" "rk4" "rk23" "rk23bs" "rk34f" [7] "rk45f" "rk45ck" "rk45e" "rk45dp6" "rk45dp7" "rk78dp" [13] "rk78f" "irk3r" "irk5r" "irk4hh" "irk6kb" "irk4l" [19] "irk6l" "ode23" "ode45" This list also contains implicit Runge-Kutta's (irk..), but they are not yet optimally coded. The only well-implemented implicit Runge-Kutta is the radau method (Hairer and Wanner 2010). --------------------------------------- References: Hairer E, Wanner G (2010). "Solving Ordinary Differential Equations II: Stiff and Differential- Algebraic Problems. Second Revised Edition. Springer-Verlag, Heidelberg. Hindmarsh, AC (1983). "ODEPACK, a Systematized Collection of ODE Solvers." In R Stepleman (ed.), "Scientifc Computing, Vol. 1 of IMACS Transactions on Scientifc Computation," pp. 55-64, IMACS / North-Holland, Amsterdam. Petzold, LR (1983). "Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations." SIAM Journal on Scientifc and Statistical Computing, 4, 136-148. Soetaert, Karline; Cash, Jeff; and Mazzia, Francesca: "Solving Differential Equations in R", 2012, Springer-Verlag. ----------------------- Ted Woollett ```
 [Maxima-discuss] RFC : Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems] From: Helmut Jarausch - 2014-04-16 19:53:58 ```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, linearly-implicit 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, Helmut See the examples and the core code at http://www.igpm.rwth-aachen.de/jarausch/Maxima As usual, there are no warranties, it's research code which is hopefully useful to someone. ```
 Re: [Maxima-discuss] Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems] From: Edwin Woollett - 2014-04-17 20:39:40 ```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, linearly-implicit 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.rwth-aachen.de/jarausch/Maxima --------------------------------------------------------------------------------- This code needs a user-friendly 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,1e-6), 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 ```
 Re: [Maxima-discuss] Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems] From: Panagiotis Papasotiriou - 2014-04-20 13:53:44 Attachments: Message as HTML ```The best way to deal with stiff ODEs would be a multi-step method able to switch between non-stiff and stiff ODE solvers automatically. As already mentioned by others here, 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 "non-stiff" mode automatically, but also has built-in 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 initial-value problem. Packages like Scilab or the R language use ODEPACK internally to solve initial-value 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 non-stiff 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 Lisp-readable 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 interfaceto 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 Maxima-readable version of ODEPACK does exist. 2014-04-17 23:39 GMT+03:00 Edwin 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, linearly-implicit 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.rwth-aachen.de/jarausch/Maxima > > --------------------------------------------------------------------------------- > > This code needs a user-friendly 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,1e-6), > 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 > _______________________________________________ > Maxima-discuss mailing list > Maxima-discuss@... > https://lists.sourceforge.net/lists/listinfo/maxima-discuss > ```
 Re: [Maxima-discuss] RFC : Solver for stiff ODEs [was:feature request: stiff ode solvers for initial value problems] From: Robert Dodier - 2014-04-17 17:46:35 ```On 2014-04-16, Helmut Jarausch wrote: > See the examples and the core code at > > http://www.igpm.rwth-aachen.de/jarausch/Maxima Thanks, Helmut. I wonder how this method compares to others for the examples you tried. Can you run any of the other methods available for Maxima on the same problems? I know there is at least one Runge-Kutta implementation, maybe other methods too. That would help us see the advantages of this new code. Also, I wonder if you will release your code under terms of the GNU General Public License, as is the rest of Maxima. It helps to have a definite statement as to what you are allowing others to do with the code. best, Robert Dodier ```