In this final blog post I want to summarize the work done so far for the GSoC. First of all I want to thank Octave community for their support:
when I had some problems I have always found someone in the mailing list or in the IRC channel willing to help me.
In my public repository I pushed three final commits:
https://bitbucket.org/Francesco_Faccio/octave-gsoc/commits/all
The first commit I pushed regards the work we did in the first part of the GSoC and it is a patch for Octave which allows to include IDA module of Sundials library as a dependence of Octave.
I wrote some configuration tests, checking for sundials_ida and sundials_nvecserial headers and functions and a macro to test whether sundials_ida has been configured with double precision realtype. I included a test for klu library in suitesparse, because it is required for sundials' sparse solver.
If all these conditions are satisfied, you can use ode15i and ode15s in Octave (you can see in my previous post how to configure sundials and octave correctly and how to run benchmark problems).
The second commit regards my work on solvers ode23 and ode45 which were already in Octave.
Since there were a lot of duplicate code, we decided to refactor the solvers, exporting the check of the options into two private functions:
$[defaults, classes, attributes] = odedefaults (n, t0, tf)$
is a function which returns a struct defaults, containing the default values of all the solvers, and structs classes and attributes, containing the arguments for validateattributes.
$options = odemergeopts (useroptions, options, classes, attributes, fun_name);$
is a function which loops over default fields and, if a user provides a new value, uses validateattributer or validatestring to test the value and overwrites the default one.
odeset.m and odeget.m have been modified in such a way that the user can add new special fields to an ODE option structure. They don't depend on private functions ode_struct_value_check and known_option_names, but use inputParser to perform input checking. This version of odeset doesn't include the option to catch a misspelled option name yet, so that 'Rel' doesn't become automatically 'RelTol'.
We wrote function decic, which can compute consistent initial conditions for ode15i, close to initial guesses y0 and yp0.
The third commit contains files ode15i.m, ode15s.m and a common interface __ode15__.cc
ode15i
The main functionalities of ode15i have been achieved. Here's a list of the options accepted by the solver:
RelTol: You can specify a relative error tolerance as a positive scalar. The default value is 1e-3.
AbsTol: You can specify an absolute error tolerance as a positive scalar or vector. The default value is 1e-6.
NormControl: Sundials controls the relative error using the norm of the solution at each step, so option NormControl is always "on".
OutputFcn, OutputSel: You can specify them to plot the components of the solution you want after each successful integration step. The default output function is $@odeplot$, but you can provide a custom output function.
Refine: You can increase the number of output points using the refine factor. __ode15__.cc uses the IDAGetDky function to reach the highest possible order of interpolation.
Stats: Only standard statistics are printed (number of successful steps, number of failed attempts and the number of function evaluations).
Events: If you provide an event function, it must be of the form $ [value,isterminal,direction] = EventFcn(t,y,yp)$. The solver localizes a change of sign on your event function and uses linear interpolation to extimate the time when the event occurs. A way to improve this functionality will be to use the rootfinding function of Sundials in order to have an higher order of interpolation.
MaxStep: Use this option to specify an upper bound for the integration step. The default value is $0.1*abs(tf - t0)$.
InitialStep: Use this option to specify the first integration step of the solver. Its default value depends on the tspan and the initial slope.
Jacobian: You can specify the Jacobian as a function which returns the modified Jacobian $[\frac{\partial F}{\partial y}, \frac{\partial F}{\partial y'}] = jac (t, y, y')$ or as a cell array containing the two constant Jacobian matrices $ \{[\frac{\partial F}{\partial y}], [\frac{\partial F}{\partial y'}] \}$.
In both cases, if the two parts of the Jacobian are sparse, the Jacobian will be stored in CCS format.
If you don't provide the Jacobian, it will be extimated internally by Sundials, using difference-quotient approximation.
In this implementation you can't provide only one part of the Jacobian yet, but you have to specify both.
JPattern, Vectorized: These options haven't been developed yet.
MaxOrder: You can specify an upper bound for the order of the solver. The default value is 5.
This solver can be extended in order to increase the compatibility with Matlab:
- completing the missing options
- allowing the user to solve equations with complex number
- allowing the user to use single precision variables.
ode15s
This solver uses the same dld-function of ode15i. Inside the m-file, the semi-explicit DAE is written as an implicit one and solve as ode15i.
The solver uses the same options of ode15i, but in addition, it has:
NonNegative: this option has not been implemented yet. It can be implemented using IDASetConstraint function.
Events: Same as ode15i, but it must have the form $[value,isterminal,direction] = EventsFcn(t,y)$.
Jacobian: You can provide a function which returns the Jacobian of the form $[\frac{\partial F}{\partial y}] = jac(t,y)$ or a constant Jacobian matrix.
Mass: You can provide a constant Mass matrix or a function handle of the form $M(t,y)$ or $M(t)$. If the Mass matrix depends on the state, the Jacobian Matrix will be computed internally from Sundials. If you provide a sparse Jacobian and a sparse Mass matrix, they will be stored in CCS format.
MStateDependence: Specify MStateDependence to "strong" or "weak" if the Mass matrix has the form $M(t,y)$. Use the value "none" if the Mass matrix has the form $M(t)$ or it is constant.
MvPattern: This option has not been developed yet.
MassSingular: This option has not been developed yet. The solver solves the problem both with and without a singular mass matrix in the same way.
InitialSlope: You can specify the initial slope of the problem. The default value is a vector of zeros.
BDF: this implementation of ode15s uses always BDF method.
This solver can be extended in order to increase the compatibility with Matlab:
- completing the missing options
- allowing the user to solve equations with complex number
- allowing the user to use single precision variables.
This code has not been merged yet, since Octave is concluding a release process.
Concluding, I really enjoied working on this project and there is still a lot of work to do. I will work in the next weeks to extend the functionalities of ode15i and ode15s. I plan to continue to contribute to Octave development, maybe both in ode and in some statistics project.
Nessun commento:
Posta un commento