Mercurial > forge
changeset 3720:45ca57d7daf5 octave-forge
Updated.
| author | treichl |
|---|---|
| date | Wed, 08 Aug 2007 20:08:00 +0000 |
| parents | 9d6873c2242b |
| children | d95f02e8540b |
| files | main/odepkg/doc/odepkg.texi |
| diffstat | 1 files changed, 638 insertions(+), 603 deletions(-) [+] |
line wrap: on
line diff
--- a/main/odepkg/doc/odepkg.texi Wed Aug 08 19:40:41 2007 +0000 +++ b/main/odepkg/doc/odepkg.texi Wed Aug 08 20:08:00 2007 +0000 @@ -1,603 +1,638 @@ -\input texinfo @c -*-texinfo-*- - -@c Copyright (c) 2007 Thomas Treichl <treichl@users.sourceforge.net> -@c odepkg.texi - A package for solving differential equations with Octave -@c -@c Permission is granted to copy, distribute and/or modify this document -@c under the terms of the GNU Free Documentation License, Version 1.2 -@c or any later version published by the Free Software Foundation; -@c with no Invariant Sections, no Front-Cover Texts, and no Back-Cover -@c Texts. A copy of the license is included in the section entitled "GNU -@c Free Documentation License". - -@c %*** Start of HEADER -@setfilename odepkg.info -@settitle OdePkg - A package for solving differential equations with Octave -@afourpaper -@c @afourwide -@c %*** End of the HEADER - -@c %*** Start of COPYRIGHT -@copying -Permission is granted to copy, distribute and/or modify this document -under the terms of the GNU Free Documentation License, Version 1.2 or -any later version published by the Free Software Foundation; with no -Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A -copy of the license is included in the section entitled "GNU Free -Documentation License". -@end copying -@c %*** End of COPYRIGHT - -@c %*** Start of TITLEPAGE -@titlepage -@title OdePkg -@subtitle A package for solving differential equations with Octave -@subtitle @b{This document currently is under development} -@author by Thomas Treichl -@page -@vskip 0pt plus 1filll -@insertcopying -@end titlepage -@c %*** End of TITLEPAGE - -@c %*** Start of BODY -@contents -@ifnottex -@node Top, Beginner's Guide, (dir), (dir) -@top Copyright -@insertcopying -@end ifnottex - -@menu -* Beginner's Guide:: Manual for users that are completely new to OdePkg -* User's Guide:: Manual for users that are already familiar with OdePkg -* Coder's Guide:: Manual for users that want to make changes to OdePkg -* Appendix:: The GNU Free Documentation License of this document -@end menu - -@c %*** Start of first chapter: Beginner's Guide -@node Beginner's Guide, User's Guide, Top, Top -@chapter Beginner's Guide -The ``Beginner's Guide'' is intended for new users who want to solve differential equations with the higher level language Octave and the toolbox OdePkg. In this chapter it will be explained what OdePkg is about in @ref{About OdePkg} and how OdePkg grew up from the beginning in @ref{OdePkg history and roadmap}. In @ref{Installation and deinstallation} it is explained how OdePkg can be installed and in @ref{First tests and demos} the first examples are explained. - -@menu -* About OdePkg:: An introduction about OdePkg -* OdePkg history and roadmap:: From the initial release until now -* Installation and deinstallation:: Setting up OdePkg on your system -* Reporting Bugs:: Writing comments and bugs -* First tests and demos:: The --foo-- example and others -@end menu - -@node About OdePkg, OdePkg history and roadmap, Beginner's Guide, Beginner's Guide -@section About OdePkg -OdePkg is part of the @b{GNU Octave Repository} (resp. the Octave--Forge project) that was initiated by Paul Kienzle in the year 2000 and that is hosted at @url{http://octave.sourceforge.net}. The package includes commands for setting up various options, output functions etc. before solving a set of differential equations with the solver functions that are also included. OdePkg formerly was initiated to solve explicitly formulated ordinary differential equations (ODEs) only, but there are already improvements so that differential algebraic equations (DAEs) in explicit form can also be solved. At this time OdePkg is under development with the main target, to make a package that is mostly compatible to commercial solver products. - -@node OdePkg history and roadmap, Installation and deinstallation, About OdePkg, Beginner's Guide -@section OdePkg history and roadmap -@cindex history -@cindex roadmap - -@multitable @columnfractions .25 .75 -@item OdePkg Version 0.0.1 -@tab The initial release was already a modification of the old ``ode package'' that was hosted at Octave--Forge and that was written by Marc Compere some when between 2000 and 2001. The four variable step--size Runge--Kutta algorithms in three solver files and the three fixed step--size solvers have been merged. It was possible to set some options for these solvers. The four output--functions (@command{odeprint}, @command{odeplot}, @command{odephas2} and @command{odephas3}) have been added along with other examples that initially have not been there. -@item OdePkg Version 0.1.x -@tab The major milestone along versions 0.1.x was that four stable solvers have been implemented (ie. @command{ode23}, @command{ode45}, @command{ode54} and @command{ode78}) supporting all options that can be set for these kind of solvers and also all necessary functions for setting their options (eg. @command{odeset}, @command{odepkg_structure_check}, etc.). Since version 0.1.3 there is also code available that interfaces the Fortran solver @file{dopri5.f} that is written by Ernst Hairer and Gerhard Wanner (cf. @file{odepkg_mexsolver_dopri5.c} and the helper files @file{odepkgext.c} and @file{odepkgmex.c}). -@item OdePkg Version 0.2.x -@tab The main work along version 0.2.x was making the interface functions for the non--stiff and stiff solvers from Ernst Hairer and Gerhard Wanner enough stable so that they could be compiled and installed by default. Wrapper functions have been added to the package with help texts and tests (eg. @command{ode2r}, @command{ode5r}, @command{oders} etc.). Six testsuite functions have been added for performance tests of the different solvers (@command{odepkg_testsuite_chemakzo}, @command{odepkg_testsuite_oregonator}, @command{odepkg_testsuite_transistor}, etc.). -@item @b{(current)} Version 0.3.x -@tab Ongoing work with this manual. Fetching and adding DDE--solvers from Netlib or Ernst Hairer or Jeff Cash if GPL compatible. -@item @b{(future)} Version 0.4.x -@tab Ongoing work with this manual. Fetching and adding IDE--solvers from Netlib or Ernst Hairer or Jeff Cash if GPL compatible. -@item @b{(future)} Version 0.5.x -@tab (Maybe) A lot of compatibility tests. -@item @b{(future)} Version 0.6.x -@tab (Maybe) Final release before version 1.0.0. -@item @b{(future)} Version 1.0.0 -@tab Completed odepkg release 1.0.0 with m--solvers and mex--solvers. -@end multitable - -@node Installation and deinstallation, Reporting Bugs, OdePkg history and roadmap, Beginner's Guide -@section Installation and deinstallation -@cindex installation -@cindex deinstallation - -OdePkg can be installed easily using the @command{pkg} command of Octave. For this get into the directory where the current release of OdePkg can be found, start @command{octave} and type -@example -pkg install odepkg-x.x.x.tar.gz -@end example -where @file{x.x.x} in the name of the @file{*.tar.gz} file is the current release number of OdePkg that is available. If you want to deinstall resp. remove OdePkg then simply type -@example -pkg uninstall odepkg-x.x.x.tar.gz -@end example -If you encounter problems during the installation process of OdePkg with the @command{pkg} command or if you have an OdePkg that seems to be broken then please report this on the mailing--list of Octave--Forge using the email address -@ifnothtml -@email{octave-dev@@lists.sourceforge.net} -@end ifnothtml -@ifhtml -@email{octave-dev @{at) lists.sourceforge.net} (replace @{at) with @@) -@end ifhtml -. - -@node Reporting Bugs, First tests and demos, Installation and deinstallation, Beginner's Guide -@section Reporting Bugs -@cindex bugs -If you encounter problems while using OdePkg or if you find bugs in the source codes then please report that via email at the Octave--Forge mailing--list using the email address -@ifnothtml -@email{octave-dev@@lists.sourceforge.net} -@end ifnothtml -@ifhtml -@email{octave-dev @{at] lists.sourceforge.net} (replace @{at] with @@) -@end ifhtml -and directly send a copy to the email address -@ifnothtml -@email{treichl@@users.sourceforge.net} -@end ifnothtml -@ifhtml -@email{treichl (at@} users.sourceforge.net} (replace (at@} with @@) -@end ifhtml -. - -@node First tests and demos, , Reporting Bugs, Beginner's Guide -@section First tests and demos -@cindex First tests and demos -Let's have a look at the first ordinary differential equation with the name ``foo''. The ``foo'' equation of second order may be of the form @math{y''(t) + C_1 y'(t) + C_2 y(t) = C_3}. With the substitutions @math{y_1(t) = y(t)} and @math{y_2(t) = y'(t)} this differential equation of second order can be split into two differential equations of first order, ie. @math{y'_1(t) = y_2(t)} and @math{y'_2(t) = - C_1 y_2(t) - C_2 y_1(t) + C_3}. Next the numerical values for the constants need to be defined, ie. @math{C_1 = 2.0}, @math{C_2 = 5.0}, @math{C_3 = 10.0}. This set of ordinary differential equations can now be written as an Octave function like -@example -function vdy = foo (vt, vy, varargin) - vdy(1,1) = vy(2); - vdy(2,1) = - 2.0 * vy(2) - 5.0 * vy(1) + 10.0; -endfunction -@end example -It is seen that this ODEs do not depend on time, nevertheless the first input argument of this function needs to be defined as the time argument, followed by a integrated state argument @command{vy} and a variable size input argument @command{varargin} that can be used to set up user defined constants or control variables. As it is known that ``foo'' is a set of @slanted{ordinary} differential equations we can choose one of the four m--file Runge--Kutta solvers (cf. @ref{Solver families}). It is also known that the time period of interest may be between @math{t_0=0.0} and @math{t_e = 5.0} as well as that the initial values of the ODEs are @math{y_1(t=0) = 0.0} and @math{y_2(t=0) = 0.0}. Solving this set of ODEs can be done by typing the following commands in the Octave interpreter window -@example -ode45 (@@foo, [0 5], [0 0]); -@end example -A figure window opens and it can be seen how this ODEs are solved from @math{t_0=0.0} to @math{t_e = 5.0}. If opening the figure window is unwanted then output arguments have to be used to catch the results of the solving process, eg. -@example -[t, y] = ode45 (@@foo, [0 5], [0 0]); -@end example -Results can also be obtained as a Octave structure if one output argument is used like in the following example. Then the results are stored in the fields @command{S.x} and @command{S.y}. -@example -S = ode45 (@@foo, [0 5], [0 0]); -@end example -As noticed before, a function for the ordinary differential equations must not be rewritten all the time if some of the parameters are going to change, that's what the input argument @command{varargin} is for. So rewrite the function @command{foo} into @command{newfoo} the following way -@example -function vdy = newfoo (vt, vy, varargin) - vdy(1,1) = vy(2); - vdy(2,1) = -varargin@{1@}*vy(2)-varargin@{2@}*vy(1)+varargin@{3@}; -endfunction -@end example -There is nothing said about the constant values anymore, but if using the following caller routine in the Octave interpreter window then the same results can be obtained as before (ie. the parameters are directly feed through from the caller routine to the function @command{newfoo}). -@example -ode45 (@@newfoo, [0 5], [0 0], 2.0, 5.0, 10.0); -@end example -The OdePkg can do much more while solving ODEs and DAEs, eg. you can set up other output functions instead of the function @command{odeplot}. So as a last example in this beginning chapter it is shown how this can be done, ie. with the command @command{odeset} -@example -A = odeset ('OutputFcn', @@odeprint); -ode45 (@@newfoo, [0 5], [0 0], A, 2.0, 5.0, 10.0); -@end example -The options structure @command{A} that can be set up with with the command @command{odeset} must always be the fourth input argument when using the ODE--solvers and the DAE--solvers (read the help files for the other solvers if there may be changes in the future). The options that can be set are described in @ref{ODE/DAE options}. - -Other examples have also been added to the OdePkg. These example files and functions are of the form @command{odepkg_equations_*}. Demos have been added to these functions as well so their functionality can be viewed very quickly. - -Before you continue reading the next chapter notice that nearly every function that belongs to the OdePkg has its own help description and demos, so you can have a look for yourself how the different functions can be used. If you want to have a look at the help description then type -@example -help fcnname -@end example -in the Octave interpreter window where @command{fcnname} is the name of the function for the help description to be viewed. Type -@example -demo fcnname -@end example -in the Octave interpreter window where @command{fcnname} is the name of the function of the demo to run. Last but not least write -@example -doc odepkg -@end example -for opening this manual in the texinfo reader of the octave interpreter window. - -@c %*** End of first chapter: Beginner's Guide - -@c %*** Start of second chapter: User's Guide -@node User's Guide, Coder's Guide, Beginner's Guide, Top -@chapter User's Guide -@cindex User's Guide -The ``User's Guide'' is intended for trained users who already do know in principal how to solve differential equations with the higher level language Octave and OdePkg. In this chapter it will be explained which solvers can be used for the different kind of problems in @ref{Solver families} and which options can be set for the optimization of the solving process in @ref{ODE/DAE options}. - -@menu -* Solver families:: The different kind of solvers -* ODE/DAE options:: Options that can be set -@end menu - -@node Solver families, ODE/DAE options, User's Guide, User's Guide -@section Solver families -@cindex Solver -In this section the different kinds of solvers are explained. It starts with the standard M--file Runge--Kutta solvers in section @ref{M--file Runge--Kutta solvers} and then it is continued with the Mex--file Hairer--Wanner solvers in section @ref{Mex--file Hairer--Wanner solvers}. Other solvers are described in section @ref{Other solvers}. Performance tests have also been added to the OdePkg. Some of these performance results have been added to section @ref{ODE solver performances}. - -@menu -* M--file Runge--Kutta solvers:: The ODE solvers written in the Octave language -* Mex--file Hairer--Wanner solvers:: Fast ODE/DAE solvers written in C--Mex -* Other solvers:: Not implemented by now -* ODE solver performances:: Cross Platform performance tests -@end menu - -@node M--file Runge--Kutta solvers, Mex--file Hairer--Wanner solvers, Solver families, Solver families -@subsection M--file Runge--Kutta solvers -The M--file Runge--Kutta solvers are written in the Octave interpreter language and are stored @file{*.m} files. There have been implemented four different solvers of similiar structure and types, ie. @command{ode23}, @command{ode45}, @command{ode54} and @command{ode78}@footnote{The descriptions for these Runge--Kutta solvers have been taken from the help texts of the initial m--file Runge--Kutta solvers that were written by Marc Compere, he also pointed out that ''a relevant discussion on step size choice can be found on page 90ff in U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 1998''.}. - -The order of all of the following Runge--Kutta methods is the order of the local truncation error, which is the principle error term in the portion of the Taylor series expansion that gets dropped, or intentionally truncated. This is different from the local error which is the difference between the estimated solution and the actual, or true solution. The local error is used in stepsize selection and may be approximated by the difference between two estimates of different order, @math{l(h) = x(O(h+1)) - x(O(h))}. With this definition, the local error will be as large as the error in the lower order method. The local truncation error is within the group of terms that gets multipled by @math{h} when solving for a solution from the general Runge--Kutta method. Therefore, the order--p solution created by the Runge--Kunge method will be roughly accurate to @math{O(h^{(p+1)})} since the local truncation error shows up in the solution as @math{e = h d}, which is @math{h} times an @math{O(h^p)} term, or rather @math{O(h^{(p+1)})}. - -@table @code -@item ode23 -Integrates a system of ordinary differential equations using second and third order Runge--Kutta formulas. This particular third order method reduces to Simpson's @math{1/3} rule and uses the third order estimation for the output solutions. Third order accurate Runge--Kutta methods have local and global errors of @math{O(h^4)} and @math{O(h^3)} respectively and yield exact results when the solution is a cubic (the variable @math{h} is the step size from one integration step to another integration step). This solver requires three function evaluations per integration step. - -@item ode45 -Integrates a system of ordinary differential equations using fourth and fifth order embedded formulas from Fehlberg. This is a fourth--order accurate integrator therefore the local error normally expected is @math{O(h^5)}. However, because this particular implementation uses the fifth--order estimate for @math{x_{out}} (ie. local extrapolation) moving forward with the fifth--order estimate should yield local error of @math{O(h^6)}. This solver requires six function evaluations per integration step. - -@item ode54 -The Fehlberg @math{4(5)} of the @command{ode45} pair is established and works well, however, the Dormand--Prince @math{4(5)} pair minimizes the local truncation error in the fifth--order estimate which is what is used to step forward (local extrapolation). Generally it produces more accurate results and costs roughly the same computationally. This solver requires seven function evaluations per integration step. - -@item ode78 -Integrates a system of ordinary differential equations using seventh and eighth order Runge--Kutta formulas. This is a seventh--order accurate integrator therefore the local error normally expected is @math{O(h^8)}. However, because this particular implementation uses the eighth--order estimate for @math{x_{out}} moving forward with the eighth--order estimate will yield errors on the order of @math{O(h^9)}. This solver requires thirteen function evaluations per integration step. - -@end table - -@node Mex--file Hairer--Wanner solvers, Other solvers, M--file Runge--Kutta solvers, Solver families -@subsection Mex--file Hairer--Wanner solvers -The mex--file Hairer--Wanner solvers are written in Fortran (hosted at @url{http://www.unige.ch/~hairer}) and have been added to the OdePkg as a compressed file with the name @file{hairer.tgz}. The licence of these solvers is a modified BSD license (without advertising clause) and can be found as @file{licence.txt} file in the @file{hairer.tgz} package and therefore the Fortran files are GPL compatible. Papers and other details about these solvers can be found at the host adress. - -Interface functions for these solvers have been created and have been added to the OdePkg. Their names are @file{odepkg_mexsolver_xxx.c} where @file{xxx} is the name of the Fortran file that is interfaced. The corresponding @file{odepkg_mexsolver_xxx.mex} files are created automatically when installing OdePkg with the @command{pkg} command, but can also be build manually with the instructions given as a preamble of every @file{odepkg_mexsolver_xxx.c} file. - -To provide a shorter name to access these solver functions also wrapper functions have been added that do link to the interface functions, eg. the command @command{oderd} links to the interface functions @command{odepkg_mexsolver_radau} and should do exactly the same. Another reason of adding wrapper functions was that help texts, demos and tests cannot be added to the @file{odepkg_mexsolver_xxx.c} files. For accessing the help texts, demos and tests for one of these solvers you should therefore always use the name of the wrapper function, eg. @command{help oderd}. - -The mex--file Hairer--Wanner solvers have been added to the OdePkg to also solve stiff ordinary differential equations that cannot be solved with one of the m--file Runge--Kutta solvers. The following table gives an overview about which solver can be used for the different kind of problems. - -@float -@multitable {ODE--Problem} {Solver name} {Wrapper file} {@file{odepkg_mexsolver_seulex.c}} {Fortran file} -@headitem ODE--Problem @tab Solver name @tab Wrapper file @tab Interface file @tab Fortran file -@item Non--stiff @tab DOPRI5 @tab @file{ode5d.m} @tab @file{odepkg_mexsolver_dopri5.c} @tab @file{dopri5.f} -@item Non--stiff @tab DOP853 @tab @file{ode8d.m} @tab @file{odepkg_mexsolver_dop853.c} @tab @file{dop853.f} -@item Non--stiff @tab ODEX @tab @file{odeox.m} @tab @file{odepkg_mexsolver_odex.c} @tab @file{odex.f} -@item Stiff @tab RADAU @tab @file{ode2r.m} @tab @file{odepkg_mexsolver_radau.c} @tab @file{radau.f} -@item Stiff @tab RADAU5 @tab @file{ode5r.m} @tab @file{odepkg_mexsolver_radau5.c} @tab @file{radau5.f} -@item Stiff @tab RODAS @tab @file{oders.m} @tab @file{odepkg_mexsolver_rodas.c} @tab @file{rodas.f} -@item Stiff @tab SEULEX @tab @file{odesx.m} @tab @file{odepkg_mexsolver_seulex.c} @tab @file{seulex.f} -@end multitable -@caption{Overview about Fortran, Interface and Wrapper files for Hairer--Wanner solvers.} -@end float - -@node Other solvers, ODE solver performances, Mex--file Hairer--Wanner solvers, Solver families -@subsection Other solvers -@c The solvers from Jeff Cash http://www.ma.ic.ac.uk/~jcash have not been added to the OdePkg because it seems that there may exist some unfixed bugs. Modified solver files with bug fixes can be found at http://pitagora.dm.uniba.it/~testset, but the disclaimer's part of this site is not clear and maybe the modified solver files are not GNU--GPL compatible. - -@c The people from the University of Bari at http://pitagora.dm.uniba.it/~testset created the solvers BIMD and GAMD. But like before the disclaimer's part of this internet site is not clear and maybe they are not GNU--GPL compatible. - -@c The OdePack solvers from http://www.netlib.org/odepack have not been added to OdePkg, because the question is not answered if they are really needed from somebody. Other solvers from http://www.netlib.org/ode have not been added because of the same reason. - -@node ODE solver performances, , Other solvers, Solver families -@subsection ODE solver performances - -@smallexample ->> odepkg ('odepkg_performance_mathires'); ------------------------------------------------------------------------------------------ - Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time ------------------------------------------------------------------------------------------ - ode113 1e-007 1e-007 1e-009 7.57 5.37 24317 21442 45760 11.697 - ode23 1e-007 1e-007 1e-009 7.23 5.03 13876 13862 41629 2.634 - ode45 1e-007 1e-007 1e-009 7.91 5.70 11017 10412 66103 2.994 - ode15s 1e-007 1e-007 1e-009 7.15 4.95 290 273 534 8 59 0.070 - ode23s 1e-007 1e-007 1e-009 6.24 4.03 702 702 2107 702 702 0.161 - ode23t 1e-007 1e-007 1e-009 6.00 3.79 892 886 1103 5 72 0.180 -ode23tb 1e-007 1e-007 1e-009 5.85 3.65 735 731 2011 5 66 0.230 ------------------------------------------------------------------------------------------ -@end smallexample - -@smallexample -octave:1> odepkg ('odepkg_performance_octavehires'); ------------------------------------------------------------------------------------------ - Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time ------------------------------------------------------------------------------------------ - ode23 1e-07 1e-07 1e-09 7.95 5.53 16179 13646 48534 168.182 - ode45 1e-07 1e-07 1e-09 8.06 5.64 9401 9398 56400 134.011 - ode54 1e-07 1e-07 1e-09 8.31 5.89 8854 7697 61971 127.261 - ode78 1e-07 1e-07 1e-09 9.06 6.64 7287 6613 94718 168.769 - odeox 1e-07 1e-07 1e-09 6.67 4.25 10969 8881 194129 226.890 - ode5d 1e-07 1e-07 1e-09 0.14 -2.28 1014 1014 6086 6.775 - ode8d 1e-07 1e-07 1e-09 0.16 -2.26 1046 1030 15385 17.602 - ode2r 1e-07 1e-07 1e-09 7.69 5.27 59 59 849 50 59 1.231 - ode5r 1e-07 1e-07 1e-09 7.55 5.13 81 81 671 71 81 1.380 - odesx 1e-07 1e-07 1e-09 6.63 4.21 39 37 1135 27 190 1.782 - oders 1e-07 1e-07 1e-09 7.08 4.66 138 138 828 138 138 2.071 ------------------------------------------------------------------------------------------ - -@end smallexample - -@smallexample ->> odepkg ('odepkg_performance_matchemakzo'); ------------------------------------------------------------------------------------------ - Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time ------------------------------------------------------------------------------------------ - ode113 1e-007 1e-007 1e-007 NaN Inf - - - - - - - ode23 1e-007 1e-007 1e-007 NaN Inf 15 15 47 0.431 - ode45 1e-007 1e-007 1e-007 NaN Inf 15 15 92 0.170 - ode15s 1e-007 1e-007 1e-007 7.04 6.20 161 154 4 35 0.521 - ode23s 1e-007 1e-007 1e-007 7.61 6.77 1676 1676 5029 1676 1677 2.704 - ode23t 1e-007 1e-007 1e-007 5.95 5.11 406 404 3 39 0.611 -ode23tb 1e-007 1e-007 1e-007 NaN Inf 607 3036 1 608 6.730 ------------------------------------------------------------------------------------------ -@end smallexample - -@smallexample -octave:1> odepkg ('odepkg_performance_octavechemakzo'); ------------------------------------------------------------------------------------------ - Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time ------------------------------------------------------------------------------------------ - ode23 1e-07 1e-07 1e-07 0.45 -0.43 432 385 1293 2.926 - ode45 1e-07 1e-07 1e-07 0.45 -0.43 277 238 1656 3.087 - ode54 1e-07 1e-07 1e-07 0.45 -0.43 216 214 1505 2.769 - ode78 1e-07 1e-07 1e-07 0.45 -0.43 210 170 2717 4.700 - ode78 1e-07 1e-07 1e-07 2.94 2.05 193 160 4815 6.150 - ode5d 1e-07 1e-07 1e-07 2.95 2.06 234 234 1406 1.499 - ode8d 1e-07 1e-07 1e-07 2.95 2.06 161 142 2056 2.149 - ode2r 1e-07 1e-07 1e-07 8.50 7.57 43 43 372 39 43 0.486 - ode5r 1e-07 1e-07 1e-07 8.50 7.57 43 43 372 39 43 0.491 - odesx 1e-07 1e-07 1e-07 7.46 6.53 22 22 502 19 96 0.597 - oders 1e-07 1e-07 1e-07 7.92 7.04 68 67 401 66 67 0.642 ------------------------------------------------------------------------------------------ -@end smallexample - -@node ODE/DAE options, , Solver families, User's Guide -@section ODE/DAE options -@cindex ODE/DAE options -The default values of odeset can be displayed if odeset is called without any input argument and one output argument argument, eg. the following way -@example -A = odeset (); -disp (A); -@end example - -@table @samp -@item RelTol -The option @option{RelTol} is used to set the relative error tolerance for the error estimation of the solver while solving. It can either be a positive scalar or a vector with every element of the vector being a positive scalar (this depends on the solver that is used). The definite error estimation equation also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T y(t), AbsTol)}. Run -@example -A = odeset ('RelTol', 1, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('RelTol', 1e-10, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{RelTol}. - -@item AbsTol -The option @option{AbsTol} is used to set the absolute error tolerance for the error estimation of the solver while solving. It can either be a positive scalar or a vector with every element of the vector being a positive scalar (this depends on the solver that is used). The definite error estimation equation also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T y(t), AbsTol)}. Run -@example -A = odeset ('AbsTol', 1e-3, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('AbsTol', 1e-10, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{AbsTol}. - -@item NormControl -The option @option{NormControl} is used to set the type of error tolerance calculation of the solver while solving. It can either be the string @command{'on'} or @command{'off'}. At the time the solver starts the initialization procedure a warning message may be displayed if the solver will ignore the @command{'on'} setting of this option because of an unhandled resp. missing implementation. The definite error estimation equation if set @command{'on'} also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T max ( norm (y(t), Inf)), AbsTol)}. Run -@example -A = odeset ('NormControl', 'on', 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('NormControl', 'off', 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{NormControl}. - -@item MaxStep -The option @option{MaxStep} is used to set the maximum step size for the solver that is used while solving. It can only be a positive scalar. By default this value is set internally by every solver and also may be different when using different solvers. Run -@example -A = odeset ('MaxStep', 10, 'OutputFcn', @@odeprint); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('MaxStep', 1e-1, 'OutputFcn', @@odeprint); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{MaxStep}. - -@item InitialStep -The option @option{InitialStep} is used to set the initial first step size for the solver. It can only be a positive scalar. By default this value is set internally by every solver and also may be different when using different solvers. Run -@example -A = odeset ('InitialStep', 1, 'OutputFcn', @@odeprint); -ode78 (@@odepkg_equations_vanderpol, [0 1], [2 0], A); -B = odeset ('InitialStep', 1e-5, 'OutputFcn', @@odeprint); -ode78 (@@odepkg_equations_vanderpol, [0 1], [2 0], B); -@end example -to see the effect of using different values for the option @option{InitialStep}. - -@item InitialSlope -The option @option{InitialSlope} is not handled by any of the solvers by now. - -@item OutputFcn -The option @option{OutputFcn} can be used to set up an output function for displaying the results of the solver while solving. It must be a function handle to a valid function. There are four predefined output functions available with OdePkg. @command{odeprint} prints the actual time values and results in the octave window while solving, @command{odeplot} plots the results over time in a new figure window while solving, @command{odephas2} plots the first result over the second result as a two--dimensional plot while solving and @command{odephas3} plots the first result over the second result over the third result as a three--dimensional plot while solving. Run -@example -A = odeset ('OutputFcn', @@odeprint); -ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); -@end example -to see the effect of using an output function with the option @option{OutputFcn}. User defined output functions can also be used. A typical framework for a self--made output function may then be of the form -@example -function [vret] = odeoutput (vt, vy, vdeci, varargin) - switch vdeci - case 'init' - ## Do everything needed to intialize output function - case 'calc' - ## Do everything needed to create output - case 'done' - ## Do everything needed to clean up output function - endswitch -endfunction -@end example -The output function @command{odeplot} is also set automatically if the solver calculation routine is called without any output argument. Run -@example -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0]); -@end example -to see an example. - -@item Refine -The option @option{Refine} is used to set the interpolation factor that is used to increase the quality for the output values if an output function is also set with the option @option{OutputFcn}. It can only be a integer value @math{0<=Refine<=5}. Run -@example -A = odeset ('Refine', 0, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('Refine', 3, 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{Refine}. - -@item OutputSel -The option @option{OutputSel} is used to set the components for which output has to be performed if an output function is also set with the option @option{OutputFcn}. It can only be a vector of integer values. Run -@example -A = odeset ('OutputSel', [1, 2], 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); -B = odeset ('OutputSel', [2], 'OutputFcn', @@odeplot); -ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); -@end example -to see the effect of using different values for the option @option{OutputSel}. - -@item Stats -The option @option{Stats} is used to print cost statistics about the solving process after solving has been finished. It can either be the string @command{'on'} or @command{'off'}. Run -@example -A = odeset ('Stats', 'off'); -[a, b] = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); -B = odeset ('Stats', 'on'); -[c, d] = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], B); -@end example -to see the effect of using different values for the option @option{Stats}. The cost statistics can also be obtained if the solver calculation routine is called with one output argument. The cost statistics then are in the output structure in field @option{stats}. Run -@example -A = odeset ('Stats', 'on'); -B = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); -disp (B); -@end example -to see an example. - -@item Jacobian -The option @option{Jacobian} can be used to set up an external Jacobian function or Jacobian matrix for DAE solvers to achieve faster and better results (ODE Runge--Kutta solvers do not need to handle a Jacobian function handle or Jacobian matrix). It must either be a function handle to a valid function or a full constant matrix of size squared the dimension of the set of differential equations. Run -@example -function vdy = fpol (vt, vy, varargin) - vdy = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; -endfunction - -function vr = fjac (vt, vy, varargin) - vr = [0, 1; ... - -1-2*vy(1)*vy(2), 1-vy(1)^2]; -endfunction - -A = odeset ('Stats', 'on'); -B = odepkg_mexsolver_radau (@@fpol, [0 20], [2 0], A); -C = odeset ('Jacobian', @@fjac, 'Stats', 'on'); -D = odepkg_mexsolver_radau (@@fpol, [0 20], [2 0], C); -@end example -to see the effect of using an Jacobian function with the option @option{Jacobian}. User defined Jacobian functions must have the form as described before (ie. @samp{function vr = fjac (vt, vy, varargin)}. - -@item JPattern -The option @option{JPattern} is not handled by any of the solvers by now. - -@item Vectorized -The option @option{Vectorized} is not handled by any of the solvers by now. - -@item Mass -The option @option{Mass} can be used to set up an external Mass function or Mass matrix for solving DAE problems. It depends on the solver that is used if @option{Mass} is supported or not. It must either be a function handle to a valid function or a full constant matrix of size squared the dimension of the set of differential equations. Run -@example -function vdy = frob (t, y, varargin) - vdy(1,1) = -0.04*y(1)+1e4*y(2)*y(3); - vdy(2,1) = 0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)^2; - vdy(3,1) = y(1)+y(2)+y(3)-1; -endfunction - -function vmas = fmas (vt, vy, varargin) - vmas = [1, 0, 0; 0, 1, 0; 0, 0, 0]; -endfunction - -A = odeset ('Mass', @@fmas); -B = oderd (@@frob, [0 1e8], [1 0 0], A); -@end example -to see the effect of using a Mass function with the option @option{Mass}. User defined Mass functions must have the form as described before (ie. @samp{function vmas = fmas (vt, vy, varargin)}. - -@item MStateDependence -The option @option{MStateDependence} can be used to set up the type of the external Mass function for solving DAE problems if a Mass function handle is set with the option @option{Mass}. It depends on the solver that is used if @option{MStateDependence} is supported or not. It must be a string of the form @command{'none'}, @command{'weak'} or @command{'strong'}. Run -@example -function vdy = frob (vt, vy, varargin) - vdy(1,1) = -0.04*vy(1)+1e4*vy(2)*vy(3); - vdy(2,1) = 0.04*vy(1)-1e4*vy(2)*vy(3)-3e7*vy(2)^2; - vdy(3,1) = vy(1)+vy(2)+vy(3)-1; -endfunction - -function vmas = fmas (vt, varargin) - vmas = [1, 0, 0; 0, 1, 0; 0, 0, 0]; -endfunction - -A = odeset ('Mass', @@fmas, 'MStateDependence', 'none'); -B = oderd (@@frob, [0 1e8], [1 0 0], A); -@end example -to see the effect of using a Mass function with the option @option{MStateDependence}. User defined Mass functions must have the form as described before (ie. @samp{function vmas = fmas (vt, varargin)} if the option @option{MStateDependence} was set to @command{'none'}, otherwise the user defined Mass function must have the form @samp{function vmas = fmas (vt, vy, varargin)} if the option @option{MStateDependence} was set to either @command{'weak'} or @command{'strong'}. - -@item MvPattern -The option @option{MvPattern} is not handled by any of the solvers by now. - -@item MassSingular -The option @option{MassSingular} is not handled by any of the solvers by now. - -@item NonNegative -The option @option{NonNegative} is not handled by any of the solvers by now. - -@item Events -The option @option{Events} can be used to set up an Event function, ie. the Event function can be used to find zero crossings in one of the results. It must either be a function handle to a valid function. Run -@example -function vdy = fbal (vt, vy, varargin) - vdy(1,1) = vy(2)+3; - vdy(2,1) = -9.81; %# m/sē -endfunction - -function [veve, vterm, vdir] = feve (vt, vy, varargin) - veve = vy(1); %# Which event component should be tread - vterm = 1; %# Terminate if an event is found - vdir = -1; %# In which direction, -1 for falling -endfunction - -A = odeset ('Events', @@feve); -B = ode78 (@@fbal, [0 1.5], [1 3], A); -plot (B.x, B.y(:,1)); -@end example -to see the effect of using an Events function with the option @option{Events}. - -@item MaxOrder -The option @option{MaxOrder} is not handled by any of the solvers by now. - -@item BDF -The option @option{BDF} is not handled by any of the solvers by now. -@end table - -@c %*** End of second chapter: User's Guide - -@c %*** Start of third chapter: Coder's Guide -@node Coder's Guide, Appendix, User's Guide, Top -@chapter Coder's Guide -@cindex Coder's Guide - -@section C Mex Function Reference -@cindex C Mex Function Reference -@include odepkgmex.texi - -@c %*** End of second chapter: Coder's Guide - -@node Appendix, , Coder's Guide, Top -@appendix -@cindex Appendix - -@menu -* GNU Free Documentation License:: -@end menu - -@include fdl.texi - -@c @node Index, C Mex Function Index, Coder's Guide, Top -@c @unnumbered Index -@c @printindex cp - -@c @node C Mex Function Index, , Coder's Guide, Top -@c @unnumbered C Mex Function Index -@c @printindex fn - -@bye - -@c Local Variables: -@c mode: texinfo -@c TeX-master: "odepkg" -@c End: +\input texinfo @c -*-texinfo-*- + +@c Copyright (c) 2007 Thomas Treichl <treichl@users.sourceforge.net> +@c odepkg.texi - A package for solving differential equations with Octave +@c +@c Permission is granted to copy, distribute and/or modify this document +@c under the terms of the GNU Free Documentation License, Version 1.2 +@c or any later version published by the Free Software Foundation; +@c with no Invariant Sections, no Front-Cover Texts, and no Back-Cover +@c Texts. A copy of the license is included in the section entitled "GNU +@c Free Documentation License". + +@c %*** Start of HEADER +@setfilename odepkg.info +@settitle OdePkg - A package for solving differential equations with Octave +@afourpaper +@c @afourwide +@c %*** End of the HEADER + +@c %*** Start of COPYRIGHT +@copying +Permission is granted to copy, distribute and/or modify this document +under the terms of the GNU Free Documentation License, Version 1.2 or +any later version published by the Free Software Foundation; with no +Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A +copy of the license is included in the section entitled "GNU Free +Documentation License". +@end copying +@c %*** End of COPYRIGHT + +@c %*** Start of TITLEPAGE +@titlepage +@title OdePkg +@subtitle A package for solving differential equations with Octave +@subtitle @b{This document currently is under development} +@author by Thomas Treichl +@page +@vskip 0pt plus 1filll +@insertcopying +@end titlepage +@c %*** End of TITLEPAGE + +@c %*** Start of BODY +@contents +@ifnottex +@node Top, Beginner's Guide, (dir), (dir) +@top Copyright +@insertcopying +@end ifnottex + +@menu +* Beginner's Guide:: Manual for users that are completely new to OdePkg +* User's Guide:: Manual for users that are already familiar with OdePkg +* Coder's Guide:: Manual for users that want to make changes to OdePkg +* Appendix:: The GNU Free Documentation License of this document +@end menu + +@c %*** Start of first chapter: Beginner's Guide +@node Beginner's Guide, User's Guide, Top, Top +@chapter Beginner's Guide +The ``Beginner's Guide'' is intended for new users who want to solve differential equations with the higher level language Octave and the package OdePkg. In this section it will be explained what OdePkg is about in @ref{About OdePkg} and how OdePkg grew up from the beginning in @ref{OdePkg history and roadmap}. In the section @ref{Installation and deinstallation} it is explained how OdePkg can be installed and in @ref{First tests and demos} the first examples are explained. + +@menu +* About OdePkg:: An introduction about OdePkg +* OdePkg history and roadmap:: From the initial release until now +* Installation and deinstallation:: Setting up OdePkg on your system +* Reporting Bugs:: Writing comments and bugs +* First tests and demos:: The --foo-- example and others +@end menu + +@node About OdePkg, OdePkg history and roadmap, Beginner's Guide, Beginner's Guide +@section About OdePkg +OdePkg is part of the @b{GNU Octave Repository} (resp. the Octave--Forge project) that was initiated by Paul Kienzle in the year 2000 and that is hosted at @url{http://octave.sourceforge.net}. The package includes commands for setting up various options, output functions etc. before solving a set of differential equations with the solver functions that are also included. OdePkg formerly was initiated to solve explicitly formulated ordinary differential equations (ODEs) only, but there are already improvements so that differential algebraic equations (DAEs) in explicit form and in implicit form (IDEs) can also be solved. At this time OdePkg is under development with the main target, to make a package that is mostly compatible to commercial solver products. + +@node OdePkg history and roadmap, Installation and deinstallation, About OdePkg, Beginner's Guide +@section OdePkg history and roadmap +@cindex history +@cindex roadmap + +@multitable @columnfractions .25 .75 +@item OdePkg Version 0.0.1 +@tab The initial release was already a modification of the old ``ode package'' that was hosted at Octave--Forge and that was written by Marc Compere some when between 2000 and 2001. The four variable step--size Runge--Kutta algorithms in three solver files and the three fixed step--size solvers have been merged. It was possible to set some options for these solvers. The four output--functions (@command{odeprint}, @command{odeplot}, @command{odephas2} and @command{odephas3}) have been added along with other examples that initially have not been there. +@item OdePkg Version 0.1.x +@tab The major milestone along versions 0.1.x was that four stable solvers have been implemented (ie. @command{ode23}, @command{ode45}, @command{ode54} and @command{ode78}) supporting all options that can be set for these kind of solvers and also all necessary functions for setting their options (eg. @command{odeset}, @command{odepkg_structure_check}, etc.). Since version 0.1.3 there is also code available that interfaces the Fortran solver @file{dopri5.f} that is written by Ernst Hairer and Gerhard Wanner (cf. @file{odepkg_mexsolver_dopri5.c} and the helper files @file{odepkgext.c} and @file{odepkgmex.c}). +@item OdePkg Version 0.2.x +@tab The main work along version 0.2.x was making the interface functions for the non--stiff and stiff solvers from Ernst Hairer and Gerhard Wanner enough stable so that they could be compiled and installed by default. Wrapper functions have been added to the package with help texts and tests (eg. @command{ode2r}, @command{ode5r}, @command{oders} etc.). Six testsuite functions have been added for performance tests of the different solvers (@command{odepkg_testsuite_chemakzo}, @command{odepkg_testsuite_oregonator}, @command{odepkg_testsuite_transistor}, etc.). +@item @b{(current)} Version 0.3.x +@tab Ongoing work with this manual. Jeff Cash released his solvers under the GPL -- first tests are done to include these solvers within OdePkg. With the beginning of version 0.3.2 new interface functions are created based on Octave's C++DLD interface to achieve a more higher performance. The first IDE--solver @file{mebdfi.f} appears and that is interfaced by @file{odepkg_octsolver.mebdfi.cc}. +@item @b{(future)} Version 0.4.x +@tab Ongoing work with this manual. Fetching and adding the DASRT IDE--solver from Netlib. Porting the mex--file solvers to Octave's C++DLD interface. +@item @b{(future)} Version 0.5.x +@tab (Maybe) A lot of compatibility tests. +@item @b{(future)} Version 0.6.x +@tab (Maybe) Final release before version 1.0.0. +@item @b{(future)} Version 1.0.0 +@tab Completed odepkg release 1.0.0 with m--solvers and DLD--solvers. +@end multitable + +@node Installation and deinstallation, Reporting Bugs, OdePkg history and roadmap, Beginner's Guide +@section Installation and deinstallation +@cindex installation +@cindex deinstallation + +OdePkg can be installed easily using the @command{pkg} command of Octave. For this get into the directory where the current release of OdePkg can be found, start @command{octave} and type +@example +pkg install odepkg-x.x.x.tar.gz +@end example +where @file{x.x.x} in the name of the @file{*.tar.gz} file is the current release number of OdePkg that is available. If you want to deinstall resp. remove OdePkg then simply type +@example +pkg uninstall odepkg-x.x.x.tar.gz +@end example +If you encounter problems during the installation process of OdePkg with the @command{pkg} command or if you have an OdePkg that seems to be broken then please report this on the mailing--list of Octave--Forge using the email address +@ifnothtml +@email{octave-dev@@lists.sourceforge.net} +@end ifnothtml +@ifhtml +@email{octave-dev @{at) lists.sourceforge.net} (replace @{at) with @@) +@end ifhtml +. + +@node Reporting Bugs, First tests and demos, Installation and deinstallation, Beginner's Guide +@section Reporting Bugs +@cindex bugs +If you encounter problems while using OdePkg or if you find bugs in the source codes then please report that via email at the Octave--Forge mailing--list using the email address +@ifnothtml +@email{octave-dev@@lists.sourceforge.net} +@end ifnothtml +@ifhtml +@email{octave-dev @{at] lists.sourceforge.net} (replace @{at] with @@) +@end ifhtml +and directly send a copy to the email address +@ifnothtml +@email{treichl@@users.sourceforge.net} +@end ifnothtml +@ifhtml +@email{treichl (at@} users.sourceforge.net} (replace (at@} with @@) +@end ifhtml +. + +@node First tests and demos, , Reporting Bugs, Beginner's Guide +@section First tests and demos +@cindex First tests and demos +Have a look at the first ordinary differential equation with the name ``foo''. The ``foo'' equation of second order may be of the form @math{y''(t) + C_1 y'(t) + C_2 y(t) = C_3}. With the substitutions @math{y_1(t) = y(t)} and @math{y_2(t) = y'(t)} this differential equation of second order can be split into two differential equations of first order, ie. @math{y'_1(t) = y_2(t)} and @math{y'_2(t) = - C_1 y_2(t) - C_2 y_1(t) + C_3}.@* +Next the numerical values for the constants need to be defined, ie. @math{C_1 = 2.0}, @math{C_2 = 5.0}, @math{C_3 = 10.0}. This set of ordinary differential equations can now be written as an Octave function like +@example +function vdy = foo (vt, vy, varargin) + vdy(1,1) = vy(2); + vdy(2,1) = - 2.0 * vy(2) - 5.0 * vy(1) + 10.0; +endfunction +@end example +It is seen that this ODEs do not depend on time nevertheless the first input argument of this function needs to be defined as the time argument followed by a integrated state argument @command{vy} as the second input argument and a variable size input argument @command{varargin} that can be used to set up user defined constants or control variables.@* +As it is known that ``foo'' is a set of @slanted{ordinary} differential equations we can choose one of the four m--file Runge--Kutta solvers (cf. @ref{Solver families}). It is also known that the time period of interest may be between @math{t_0=0.0} and @math{t_e = 5.0} as well as that the initial values of the ODEs are @math{y_1(t=0) = 0.0} and @math{y_2(t=0) = 0.0}. Solving this set of ODEs can be done by typing the following commands in the Octave interpreter window +@example +ode45 (@@foo, [0 5], [0 0]); +@end example +A figure window opens and it can be seen how this ODEs are solved from @math{t_0=0.0} to @math{t_e = 5.0}. If opening the figure window is unwanted then output arguments have to be used to catch the results of the solving process and to not pass the results to the window plotter, eg. +@example +[t, y] = ode45 (@@foo, [0 5], [0 0]); +@end example +Results can also be obtained as an Octave structure if one output argument is used like in the following example. Then the results are stored in the fields @command{S.x} and @command{S.y}. +@example +S = ode45 (@@foo, [0 5], [0 0]); +@end example +As noticed before, a function for the ordinary differential equations must not be rewritten all the time if some of the parameters are going to change, that's what the input argument @command{varargin} is for. So rewrite the function @command{foo} into @command{newfoo} the following way +@example +function vdy = newfoo (vt, vy, varargin) + vdy(1,1) = vy(2); + vdy(2,1) = -varargin@{1@}*vy(2)-varargin@{2@}*vy(1)+varargin@{3@}; +endfunction +@end example +There is nothing said about the constant values anymore, but if using the following caller routine in the Octave interpreter window then the same results can be obtained with the new function @command{newfoo} as before with the function @command{foo} (ie. the parameters are directly feed through from the caller routine @command{ode45} to the function @command{newfoo}). +@example +ode45 (@@newfoo, [0 5], [0 0], 2.0, 5.0, 10.0); +@end example +The OdePkg can do much more while solving ODEs and DAEs, eg. setting up other output functions instead of the function @command{odeplot}. So as a last example in this beginning chapter it is shown how this can be done, ie. with the command @command{odeset} +@example +A = odeset ('OutputFcn', @@odeprint); +ode45 (@@newfoo, [0 5], [0 0], A, 2.0, 5.0, 10.0); +@end example +The options structure @command{A} that can be set up with with the command @command{odeset} must always be the fourth input argument when using the ODE--solvers and the DAE--solvers but if you are using an IDE--solver then @command{A} must be the fifth input argument (read the help files for the other solvers if there may be changes in the future). The options that can be set are described in @ref{ODE/DAE/IDE options}. + +Further examples have also been added to the OdePkg. These example files and functions are of the form @command{odepkg_equations_*}. Different testsuite examples have been added to OdePkg that are stored in files with filenames @command{odepkg_testsuite_*}. Before continuing reading the next chapter note that nearly every function that comes with OdePkg has its own help description and demos. Look for yourself how the different functions, options and combinations can be used. If you want to have a look at the help description then type +@example +help fcnname +@end example +in the Octave command window where @command{fcnname} is the name of the function for the help description to be viewed. Type +@example +demo fcnname +@end example +in the Octave command window where @command{fcnname} is the name of the function of the demo to run. Last but not least write +@example +doc odepkg +@end example +for opening this manual in the texinfo reader of the octave command window. + +@c %*** End of first chapter: Beginner's Guide + +@c %*** Start of second chapter: User's Guide +@node User's Guide, Coder's Guide, Beginner's Guide, Top +@chapter User's Guide +@cindex User's Guide +The ``User's Guide'' is intended for trained users who already do know in principal how to solve differential equations with the higher level language Octave and OdePkg. In this chapter it will be explained which solvers can be used for the different kind of problems in @ref{Solver families} and which options can be set for the optimization of the solving process in @ref{ODE/DAE/IDE options}. + +@menu +* Differential Equation Problems:: Problems be solved with OdePkg +* Solver families:: Different kind of solvers +* ODE/DAE/IDE options:: Options that can be set +@end menu + +@node Differential Equation Problems, Solver families, User's Guide, User's Guide +@section Differential Equation Problems +@cindex Problems +In this section the different kind of differential equation problems are explained that can be solved with OdePkg. The formulation of ordinary differential equations is described in section @ref{ODE problems} followed by the description of explicetly formulated differential algebraic equations in section @ref{DAE problems} and implicetely formulated differential algebraic equations @ref{IDE problems}. + +@menu +* ODE problems:: Ordinary differential equations +* DAE problems:: Explicitely formulated differential algebraic equations +* IDE problems:: Implecitely formulated differential algebraic equations +@end menu + +@node ODE problems, DAE problems, Differential Equation Problems, Differential Equation Problems +@subsection ODE problems +ODE problems in general are of the form @math{y'(t) = f(t,y)} where @math{y'(t)} may be a scalar or vector of state variables. The variable @math{t} always is a scalar describing one point of time and the variable @math{y(t)} is a scalar or vector of solutions from the set of ordinary differential equations. + +@node DAE problems, IDE problems, ODE problems, Differential Equation Problems +@subsection DAE problems +DAE problems in general are of the form @math{M(t,y) \cdot y'(t) = f(t,y)} where @math{y'(t)} may be a scalar or vector of state variables. The variable @math{t} always is a scalar describing one point of time and the variable @math{y(t)} is a scalar or vector of solutions from the set of ordinary differential equations. The variable @math{M(t,y)} is the squared mass matrix that may depend on @math{y} and @math{t}. + +@node IDE problems, , DAE problems, Differential Equation Problems +@subsection IDE problems + +@node Solver families, ODE/DAE/IDE options, Differential Equation Problems, User's Guide +@section Solver families +@cindex Solver +In this section different kind of solvers are explained that have been implemented in OdePkg. This section starts with the standard M--file Runge--Kutta solvers in section @ref{M--file Runge--Kutta solvers} and is continued with the Mex--file Hairer--Wanner solvers in section @ref{Mex--file Hairer--Wanner solvers}. @c Other solvers are described in section @ref{Other solvers}. +Performance tests have also been added to the OdePkg. Some of these performance results have been added to section @ref{ODE solver performances}. + +@menu +* M--file Runge--Kutta solvers:: The ODE solvers written in the Octave language +* Mex--file Hairer--Wanner solvers:: Fast ODE/DAE solvers written in C--Mex +* Other solvers:: Not implemented by now +* ODE solver performances:: Cross Platform performance tests +@end menu + +@node M--file Runge--Kutta solvers, Mex--file Hairer--Wanner solvers, Solver families, Solver families +@subsection M--file Runge--Kutta solvers +The M--file Runge--Kutta solvers are written in the Octave interpreter language and are stored @file{*.m} files. There have been implemented four different solvers of similiar structure and types, ie. @command{ode23}, @command{ode45}, @command{ode54} and @command{ode78}@footnote{The descriptions for these Runge--Kutta solvers have been taken from the help texts of the initial m--file Runge--Kutta solvers that were written by Marc Compere, he also pointed out that ''a relevant discussion on step size choice can be found on page 90ff in U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 1998''.}. + +The order of all of the following Runge--Kutta methods is the order of the local truncation error, which is the principle error term in the portion of the Taylor series expansion that gets dropped, or intentionally truncated. This is different from the local error which is the difference between the estimated solution and the actual, or true solution. The local error is used in stepsize selection and may be approximated by the difference between two estimates of different order, @math{l(h) = x(O(h+1)) - x(O(h))}. With this definition, the local error will be as large as the error in the lower order method. The local truncation error is within the group of terms that gets multipled by @math{h} when solving for a solution from the general Runge--Kutta method. Therefore, the order--p solution created by the Runge--Kunge method will be roughly accurate to @math{O(h^{(p+1)})} since the local truncation error shows up in the solution as @math{e = h d}, which is @math{h} times an @math{O(h^p)} term, or rather @math{O(h^{(p+1)})}. + +@table @code +@item ode23 +Integrates a system of ordinary differential equations using second and third order Runge--Kutta formulas. This particular third order method reduces to Simpson's @math{1/3} rule and uses the third order estimation for the output solutions. Third order accurate Runge--Kutta methods have local and global errors of @math{O(h^4)} and @math{O(h^3)} respectively and yield exact results when the solution is a cubic (the variable @math{h} is the step size from one integration step to another integration step). This solver requires three function evaluations per integration step. + +@item ode45 +Integrates a system of ordinary differential equations using fourth and fifth order embedded formulas from Fehlberg. This is a fourth--order accurate integrator therefore the local error normally expected is @math{O(h^5)}. However, because this particular implementation uses the fifth--order estimate for @math{x_{out}} (ie. local extrapolation) moving forward with the fifth--order estimate should yield local error of @math{O(h^6)}. This solver requires six function evaluations per integration step. + +@item ode54 +The Fehlberg @math{4(5)} of the @command{ode45} pair is established and works well, however, the Dormand--Prince @math{4(5)} pair minimizes the local truncation error in the fifth--order estimate which is what is used to step forward (local extrapolation). Generally it produces more accurate results and costs roughly the same computationally. This solver requires seven function evaluations per integration step. + +@item ode78 +Integrates a system of ordinary differential equations using seventh and eighth order Runge--Kutta formulas. This is a seventh--order accurate integrator therefore the local error normally expected is @math{O(h^8)}. However, because this particular implementation uses the eighth--order estimate for @math{x_{out}} moving forward with the eighth--order estimate will yield errors on the order of @math{O(h^9)}. This solver requires thirteen function evaluations per integration step. + +@end table + +@node Mex--file Hairer--Wanner solvers, Other solvers, M--file Runge--Kutta solvers, Solver families +@subsection Mex--file Hairer--Wanner solvers +The mex--file Hairer--Wanner solvers are written in Fortran (hosted at @url{http://www.unige.ch/~hairer}) and have been added to the OdePkg as a compressed file with the name @file{hairer.tgz}. The licence of these solvers is a modified BSD license (without advertising clause) and can be found as @file{licence.txt} file in the @file{hairer.tgz} package and therefore the Fortran files are GPL compatible. Papers and other details about these solvers can be found at the host adress. + +Interface functions for these solvers have been created and have been added to the OdePkg. Their names are @file{odepkg_mexsolver_xxx.c} where @file{xxx} is the name of the Fortran file that is interfaced. The corresponding @file{odepkg_mexsolver_xxx.mex} files are created automatically when installing OdePkg with the @command{pkg} command, but can also be build manually with the instructions given as a preamble of every @file{odepkg_mexsolver_xxx.c} file. + +To provide a shorter name to access these solver functions also wrapper functions have been added that do link to the interface functions, eg. the command @command{oderd} links to the interface functions @command{odepkg_mexsolver_radau} and should do exactly the same. Another reason of adding wrapper functions was that help texts, demos and tests cannot be added to the @file{odepkg_mexsolver_xxx.c} files. For accessing the help texts, demos and tests for one of these solvers you should therefore always use the name of the wrapper function, eg. @command{help oderd}. + +The mex--file Hairer--Wanner solvers have been added to the OdePkg to also solve stiff ordinary differential equations that cannot be solved with one of the m--file Runge--Kutta solvers. The following table gives an overview about which solver can be used for the different kind of problems. + +@float +@multitable {ODE--Problem} {Solver name} {Wrapper file} {@file{odepkg_mexsolver_seulex.c}} {Fortran file} +@headitem ODE--Problem @tab Solver name @tab Wrapper file @tab Interface file @tab Fortran file +@item Non--stiff @tab DOPRI5 @tab @file{ode5d.m} @tab @file{odepkg_mexsolver_dopri5.c} @tab @file{dopri5.f} +@item Non--stiff @tab DOP853 @tab @file{ode8d.m} @tab @file{odepkg_mexsolver_dop853.c} @tab @file{dop853.f} +@item Non--stiff @tab ODEX @tab @file{odeox.m} @tab @file{odepkg_mexsolver_odex.c} @tab @file{odex.f} +@item Stiff @tab RADAU @tab @file{ode2r.m} @tab @file{odepkg_mexsolver_radau.c} @tab @file{radau.f} +@item Stiff @tab RADAU5 @tab @file{ode5r.m} @tab @file{odepkg_mexsolver_radau5.c} @tab @file{radau5.f} +@item Stiff @tab RODAS @tab @file{oders.m} @tab @file{odepkg_mexsolver_rodas.c} @tab @file{rodas.f} +@item Stiff @tab SEULEX @tab @file{odesx.m} @tab @file{odepkg_mexsolver_seulex.c} @tab @file{seulex.f} +@end multitable +@caption{Overview about Fortran, Interface and Wrapper files for Hairer--Wanner solvers.} +@end float + +@node Other solvers, ODE solver performances, Mex--file Hairer--Wanner solvers, Solver families +@subsection Other solvers +@c The solvers from Jeff Cash http://www.ma.ic.ac.uk/~jcash have not been added to the OdePkg because it seems that there may exist some unfixed bugs. Modified solver files with bug fixes can be found at http://pitagora.dm.uniba.it/~testset, but the disclaimer's part of this site is not clear and maybe the modified solver files are not GNU--GPL compatible. + +@c The people from the University of Bari at http://pitagora.dm.uniba.it/~testset created the solvers BIMD and GAMD. But like before the disclaimer's part of this internet site is not clear and maybe they are not GNU--GPL compatible. + +@c The OdePack solvers from http://www.netlib.org/odepack have not been added to OdePkg, because the question is not answered if they are really needed from somebody. Other solvers from http://www.netlib.org/ode have not been added because of the same reason. + +@node ODE solver performances, , Other solvers, Solver families +@subsection ODE solver performances + +@smallexample +>> odepkg ('odepkg_performance_mathires'); +----------------------------------------------------------------------------------------- + Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time +----------------------------------------------------------------------------------------- + ode113 1e-007 1e-007 1e-009 7.57 5.37 24317 21442 45760 11.697 + ode23 1e-007 1e-007 1e-009 7.23 5.03 13876 13862 41629 2.634 + ode45 1e-007 1e-007 1e-009 7.91 5.70 11017 10412 66103 2.994 + ode15s 1e-007 1e-007 1e-009 7.15 4.95 290 273 534 8 59 0.070 + ode23s 1e-007 1e-007 1e-009 6.24 4.03 702 702 2107 702 702 0.161 + ode23t 1e-007 1e-007 1e-009 6.00 3.79 892 886 1103 5 72 0.180 +ode23tb 1e-007 1e-007 1e-009 5.85 3.65 735 731 2011 5 66 0.230 +----------------------------------------------------------------------------------------- +@end smallexample + +@smallexample +octave:1> odepkg ('odepkg_performance_octavehires'); +----------------------------------------------------------------------------------------- + Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time +----------------------------------------------------------------------------------------- + ode23 1e-07 1e-07 1e-09 7.95 5.53 16179 13646 48534 168.182 + ode45 1e-07 1e-07 1e-09 8.06 5.64 9401 9398 56400 134.011 + ode54 1e-07 1e-07 1e-09 8.31 5.89 8854 7697 61971 127.261 + ode78 1e-07 1e-07 1e-09 9.06 6.64 7287 6613 94718 168.769 + odeox 1e-07 1e-07 1e-09 6.67 4.25 10969 8881 194129 226.890 + ode5d 1e-07 1e-07 1e-09 0.14 -2.28 1014 1014 6086 6.775 + ode8d 1e-07 1e-07 1e-09 0.16 -2.26 1046 1030 15385 17.602 + ode2r 1e-07 1e-07 1e-09 7.69 5.27 59 59 849 50 59 1.231 + ode5r 1e-07 1e-07 1e-09 7.55 5.13 81 81 671 71 81 1.380 + odesx 1e-07 1e-07 1e-09 6.63 4.21 39 37 1135 27 190 1.782 + oders 1e-07 1e-07 1e-09 7.08 4.66 138 138 828 138 138 2.071 +----------------------------------------------------------------------------------------- + +@end smallexample + +@smallexample +>> odepkg ('odepkg_performance_matchemakzo'); +----------------------------------------------------------------------------------------- + Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time +----------------------------------------------------------------------------------------- + ode113 1e-007 1e-007 1e-007 NaN Inf - - - - - - + ode23 1e-007 1e-007 1e-007 NaN Inf 15 15 47 0.431 + ode45 1e-007 1e-007 1e-007 NaN Inf 15 15 92 0.170 + ode15s 1e-007 1e-007 1e-007 7.04 6.20 161 154 4 35 0.521 + ode23s 1e-007 1e-007 1e-007 7.61 6.77 1676 1676 5029 1676 1677 2.704 + ode23t 1e-007 1e-007 1e-007 5.95 5.11 406 404 3 39 0.611 +ode23tb 1e-007 1e-007 1e-007 NaN Inf 607 3036 1 608 6.730 +----------------------------------------------------------------------------------------- +@end smallexample + +@smallexample +octave:1> odepkg ('odepkg_performance_octavechemakzo'); +----------------------------------------------------------------------------------------- + Solver RelTol AbsTol Init Mescd Scd Steps Accept FEval JEval LUdec Time +----------------------------------------------------------------------------------------- + ode23 1e-07 1e-07 1e-07 0.45 -0.43 432 385 1293 2.926 + ode45 1e-07 1e-07 1e-07 0.45 -0.43 277 238 1656 3.087 + ode54 1e-07 1e-07 1e-07 0.45 -0.43 216 214 1505 2.769 + ode78 1e-07 1e-07 1e-07 0.45 -0.43 210 170 2717 4.700 + ode78 1e-07 1e-07 1e-07 2.94 2.05 193 160 4815 6.150 + ode5d 1e-07 1e-07 1e-07 2.95 2.06 234 234 1406 1.499 + ode8d 1e-07 1e-07 1e-07 2.95 2.06 161 142 2056 2.149 + ode2r 1e-07 1e-07 1e-07 8.50 7.57 43 43 372 39 43 0.486 + ode5r 1e-07 1e-07 1e-07 8.50 7.57 43 43 372 39 43 0.491 + odesx 1e-07 1e-07 1e-07 7.46 6.53 22 22 502 19 96 0.597 + oders 1e-07 1e-07 1e-07 7.92 7.04 68 67 401 66 67 0.642 +----------------------------------------------------------------------------------------- +@end smallexample + +@node ODE/DAE/IDE options, , Solver families, User's Guide +@section ODE/DAE/IDE options +@cindex ODE/DAE/IDE options +The default values of odeset can be displayed if odeset is called without any input argument and one output argument argument, eg. the following way +@example +A = odeset (); +disp (A); +@end example + +@table @samp +@item RelTol +The option @option{RelTol} is used to set the relative error tolerance for the error estimation of the solver while solving. It can either be a positive scalar or a vector with every element of the vector being a positive scalar (this depends on the solver that is used). The definite error estimation equation also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T y(t), AbsTol)}. Run +@example +A = odeset ('RelTol', 1, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('RelTol', 1e-10, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{RelTol}. + +@item AbsTol +The option @option{AbsTol} is used to set the absolute error tolerance for the error estimation of the solver while solving. It can either be a positive scalar or a vector with every element of the vector being a positive scalar (this depends on the solver that is used). The definite error estimation equation also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T y(t), AbsTol)}. Run +@example +A = odeset ('AbsTol', 1e-3, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('AbsTol', 1e-10, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{AbsTol}. + +@item NormControl +The option @option{NormControl} is used to set the type of error tolerance calculation of the solver while solving. It can either be the string @command{'on'} or @command{'off'}. At the time the solver starts the initialization procedure a warning message may be displayed if the solver will ignore the @command{'on'} setting of this option because of an unhandled resp. missing implementation. The definite error estimation equation if set @command{'on'} also depends on the solver that is used, but generalized it may be of the form @math{e(t) = max (RelTol^T max ( norm (y(t), Inf)), AbsTol)}. Run +@example +A = odeset ('NormControl', 'on', 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('NormControl', 'off', 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{NormControl}. + +@item MaxStep +The option @option{MaxStep} is used to set the maximum step size for the solver that is used while solving. It can only be a positive scalar. By default this value is set internally by every solver and also may be different when using different solvers. Run +@example +A = odeset ('MaxStep', 10, 'OutputFcn', @@odeprint); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('MaxStep', 1e-1, 'OutputFcn', @@odeprint); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{MaxStep}. + +@item InitialStep +The option @option{InitialStep} is used to set the initial first step size for the solver. It can only be a positive scalar. By default this value is set internally by every solver and also may be different when using different solvers. Run +@example +A = odeset ('InitialStep', 1, 'OutputFcn', @@odeprint); +ode78 (@@odepkg_equations_vanderpol, [0 1], [2 0], A); +B = odeset ('InitialStep', 1e-5, 'OutputFcn', @@odeprint); +ode78 (@@odepkg_equations_vanderpol, [0 1], [2 0], B); +@end example +to see the effect of using different values for the option @option{InitialStep}. + +@item InitialSlope +The option @option{InitialSlope} is not handled by any of the solvers by now. + +@item OutputFcn +The option @option{OutputFcn} can be used to set up an output function for displaying the results of the solver while solving. It must be a function handle to a valid function. There are four predefined output functions available with OdePkg. @command{odeprint} prints the actual time values and results in the octave window while solving, @command{odeplot} plots the results over time in a new figure window while solving, @command{odephas2} plots the first result over the second result as a two--dimensional plot while solving and @command{odephas3} plots the first result over the second result over the third result as a three--dimensional plot while solving. Run +@example +A = odeset ('OutputFcn', @@odeprint); +ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); +@end example +to see the effect of using an output function with the option @option{OutputFcn}. User defined output functions can also be used. A typical framework for a self--made output function may then be of the form +@example +function [vret] = odeoutput (vt, vy, vdeci, varargin) + switch vdeci + case 'init' + ## Do everything needed to intialize output function + case 'calc' + ## Do everything needed to create output + case 'done' + ## Do everything needed to clean up output function + endswitch +endfunction +@end example +The output function @command{odeplot} is also set automatically if the solver calculation routine is called without any output argument. Run +@example +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0]); +@end example +to see an example. + +@item Refine +The option @option{Refine} is used to set the interpolation factor that is used to increase the quality for the output values if an output function is also set with the option @option{OutputFcn}. It can only be a integer value @math{0<=Refine<=5}. Run +@example +A = odeset ('Refine', 0, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('Refine', 3, 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{Refine}. + +@item OutputSel +The option @option{OutputSel} is used to set the components for which output has to be performed if an output function is also set with the option @option{OutputFcn}. It can only be a vector of integer values. Run +@example +A = odeset ('OutputSel', [1, 2], 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], A); +B = odeset ('OutputSel', [2], 'OutputFcn', @@odeplot); +ode78 (@@odepkg_equations_vanderpol, [0 20], [2 0], B); +@end example +to see the effect of using different values for the option @option{OutputSel}. + +@item Stats +The option @option{Stats} is used to print cost statistics about the solving process after solving has been finished. It can either be the string @command{'on'} or @command{'off'}. Run +@example +A = odeset ('Stats', 'off'); +[a, b] = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); +B = odeset ('Stats', 'on'); +[c, d] = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], B); +@end example +to see the effect of using different values for the option @option{Stats}. The cost statistics can also be obtained if the solver calculation routine is called with one output argument. The cost statistics then are in the output structure in field @option{stats}. Run +@example +A = odeset ('Stats', 'on'); +B = ode78 (@@odepkg_equations_vanderpol, [0 2], [2 0], A); +disp (B); +@end example +to see an example. + +@item Jacobian +The option @option{Jacobian} can be used to set up an external Jacobian function or Jacobian matrix for DAE solvers to achieve faster and better results (ODE Runge--Kutta solvers do not need to handle a Jacobian function handle or Jacobian matrix). It must either be a function handle to a valid function or a full constant matrix of size squared the dimension of the set of differential equations. Run +@example +function vdy = fpol (vt, vy, varargin) + vdy = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; +endfunction + +function vr = fjac (vt, vy, varargin) + vr = [0, 1; ... + -1-2*vy(1)*vy(2), 1-vy(1)^2]; +endfunction + +A = odeset ('Stats', 'on'); +B = odepkg_mexsolver_radau (@@fpol, [0 20], [2 0], A); +C = odeset ('Jacobian', @@fjac, 'Stats', 'on'); +D = odepkg_mexsolver_radau (@@fpol, [0 20], [2 0], C); +@end example +to see the effect of using an Jacobian function with the option @option{Jacobian}. User defined Jacobian functions must have the form as described before (ie. @samp{function vr = fjac (vt, vy, varargin)}. + +@item JPattern +The option @option{JPattern} is not handled by any of the solvers by now. + +@item Vectorized +The option @option{Vectorized} is not handled by any of the solvers by now. + +@item Mass +The option @option{Mass} can be used to set up an external Mass function or Mass matrix for solving DAE problems. It depends on the solver that is used if @option{Mass} is supported or not. It must either be a function handle to a valid function or a full constant matrix of size squared the dimension of the set of differential equations. Run +@example +function vdy = frob (t, y, varargin) + vdy(1,1) = -0.04*y(1)+1e4*y(2)*y(3); + vdy(2,1) = 0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)^2; + vdy(3,1) = y(1)+y(2)+y(3)-1; +endfunction + +function vmas = fmas (vt, vy, varargin) + vmas = [1, 0, 0; 0, 1, 0; 0, 0, 0]; +endfunction + +A = odeset ('Mass', @@fmas); +B = oderd (@@frob, [0 1e8], [1 0 0], A); +@end example +to see the effect of using a Mass function with the option @option{Mass}. User defined Mass functions must have the form as described before (ie. @samp{function vmas = fmas (vt, vy, varargin)}. + +@item MStateDependence +The option @option{MStateDependence} can be used to set up the type of the external Mass function for solving DAE problems if a Mass function handle is set with the option @option{Mass}. It depends on the solver that is used if @option{MStateDependence} is supported or not. It must be a string of the form @command{'none'}, @command{'weak'} or @command{'strong'}. Run +@example +function vdy = frob (vt, vy, varargin) + vdy(1,1) = -0.04*vy(1)+1e4*vy(2)*vy(3); + vdy(2,1) = 0.04*vy(1)-1e4*vy(2)*vy(3)-3e7*vy(2)^2; + vdy(3,1) = vy(1)+vy(2)+vy(3)-1; +endfunction + +function vmas = fmas (vt, varargin) + vmas = [1, 0, 0; 0, 1, 0; 0, 0, 0]; +endfunction + +A = odeset ('Mass', @@fmas, 'MStateDependence', 'none'); +B = oderd (@@frob, [0 1e8], [1 0 0], A); +@end example +to see the effect of using a Mass function with the option @option{MStateDependence}. User defined Mass functions must have the form as described before (ie. @samp{function vmas = fmas (vt, varargin)} if the option @option{MStateDependence} was set to @command{'none'}, otherwise the user defined Mass function must have the form @samp{function vmas = fmas (vt, vy, varargin)} if the option @option{MStateDependence} was set to either @command{'weak'} or @command{'strong'}. + +@item MvPattern +The option @option{MvPattern} is not handled by any of the solvers by now. + +@item MassSingular +The option @option{MassSingular} is not handled by any of the solvers by now. + +@item NonNegative +The option @option{NonNegative} can be used to set single solution variables to zero even if their real solution would be a negative value. It must be a vector describing the positions in the solution vector for which the option @option{NonNegative} should be used. Run +@example +vfun = @@(vt,vy) -abs(vy); +vopt = odeset ('NonNegative', [1]); + +[vt1, vy1] = ode78 (vfun, [0 100], [1]); +[vt2, vy2] = ode78 (vfun, [0 100], [1], vopt); + +subplot (2,1,1); plot (vt1, vy1); +subplot (2,1,2); plot (vt2, vy2); +@end example +to see the effect of not using the option @option{NonNegative} in the upper subplot and if using the option @option{NonNegative} in the lower suplot. + +@item Events +The option @option{Events} can be used to set up an Event function, ie. the Event function can be used to find zero crossings in one of the results. It must either be a function handle to a valid function. Run +@example +function vdy = fbal (vt, vy, varargin) + vdy(1,1) = vy(2)+3; + vdy(2,1) = -9.81; %# m/sē +endfunction + +function [veve, vterm, vdir] = feve (vt, vy, varargin) + veve = vy(1); %# Which event component should be tread + vterm = 1; %# Terminate if an event is found + vdir = -1; %# In which direction, -1 for falling +endfunction + +A = odeset ('Events', @@feve); +B = ode78 (@@fbal, [0 1.5], [1 3], A); +plot (B.x, B.y(:,1)); +@end example +to see the effect of using an Events function with the option @option{Events}. + +@item MaxOrder +The option @option{MaxOrder} is not handled by any of the solvers by now. + +@item BDF +The option @option{BDF} is not handled by any of the solvers by now. +@end table + +@c %*** End of second chapter: User's Guide + +@c %*** Start of third chapter: Coder's Guide +@node Coder's Guide, Appendix, User's Guide, Top +@chapter Coder's Guide +@cindex Coder's Guide + +@section C Mex Function Reference +@cindex C Mex Function Reference +@include odepkgmex.texi + +@c %*** End of second chapter: Coder's Guide + +@node Appendix, , Coder's Guide, Top +@appendix +@cindex Appendix + +@menu +* GNU Free Documentation License:: +@end menu + +@include fdl.texi + +@c @node Index, C Mex Function Index, Coder's Guide, Top +@c @unnumbered Index +@c @printindex cp + +@c @node C Mex Function Index, , Coder's Guide, Top +@c @unnumbered C Mex Function Index +@c @printindex fn + +@bye + +@c Local Variables: +@c mode: texinfo +@c TeX-master: "odepkg" +@c End: