martedì 23 agosto 2016

Final Evaluation

Dear all,

GSoC16 is at the end, I want to thank you for giving me the opportunity to participate in this unique experience.

My goal was to add two solvers to the odepkg package, two exponential integrators, that are a class of numerical methods for the solution of partial and ordinary differential equations.
 
  • Here is a summary of the work done in this months:

For the midterm evaluation (repository of the midterm work): 

After a first period of study, I started with implementations.
I start with the implementation of the φ functions (phi1m.m, phi2m.m, phi3m.m, phi4m.m), necessary to calculate the matrix functions in the two schemes (exprk.m, exprb.m) for a general exponential Runge-Kutta and Rosenbrock integrator. These schemes have been very useful to verify the correctness of the official methods. I applied them to four different exponential methods (exprk2.m, exprk3.m, exprb3.m, exprb4.m).
Then I wrote a code that implements the augmented matrix  Ã (augmat.m) described in theorem 2.1 of [HAM 11]. This matrix will be given as input to expmv and will allow us to implement the exponential integrators by evaluating a single exponential of this augmented matrix avoiding the need to compute any φ functions.

For the final evaluation (REPOSITORY OF THE FINAL WORK):

After the midterm evaluation, I implemented the official exponential method that will be in odepkg, it is the advanced third order method with adaptive time stepping from Runge-Kutta family. Unlike the general schemes, for this code I used the function expmv so that I could evaluate a single exponential of an augmented matrix instead of computing any φ functions. I wrote a code for the stepper (exp_runge_kutta_23) and for the solver (odexprk23), according to odepkg standards.
Then I changed the file augmat.m to calculate a "double" augmented matrix, and in this way, I have been able to calculate, with a single call to expmv, both the solution to the next step, and the lower order solution for the estimation of the error.
In the same way I implemented the official exponential fourth order method with adaptive time stepping from Rosenbrock family (exp_rosenbrock_34, odexprb34)
I finally tested the order of my methods and I translate in C++ a for loop (__taylor__.cc) of the Nick Higham function expmv to improve the speed of an algorithm.

  • In the repository of the final work you will find a patch for odepkg.

I added to the package the following files

/inst/odexprk23.m : Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a third-order exponential Runge-Kutta method.
/inst/odexprb34.m : Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a fourth-order exponential Rosenbrock method.
/inst/steppers/exp_runge_kutta_23.m : Stepper for the odexprk23 solver.
/inst/steppers/exp_rosenbrock_34.m : Stepper for the odexprb34 solver.
/inst/utilities/expmv.m : Compute the matrix exponential times vector or matrix (2010 Awad H. Al-Mohy and Nicholas J. Higham).
/inst/utilities/augmat.m, estimates.m, expmv_tspan.m, select_taylor_degree_aug.m, NormAm.m, select_taylor_degree.m, theta_taylor, theta_taylor_half, theta_taylor_single : Functions related to expmv.
/src/__taylor__.cc : for loop of expmv. 

  • I export the patch from
odepkg changeset: 542:3f74a6e911e1,
octave version: 4.1.0+, changeset: 22348:9deb86bb5632.
My patch exponential_integrators.diff requires a changeset >= 22345.
To apply the patch follow the commands

hg import exponential_integrators.diff
hg archive ../odepkg-0.9.1.tar.gz

run Octave 4.1.0+

pkg install odepkg-0.9.1.tar.gz
pkg load odepkg 

to test the patch:


demo odexprk23
demo odexprb34 (slightly slower)

if everything goes well, the demo gives a plot on convergence of methods.

  •  Demos give the following warnings
warning: odexprk23: unknown field 'LinOp' in ODE options
warning: odexprk23: unknown field 'GFun' in ODE options

warning: odexprb34: unknown field 'DFdt' in ODE options  
I haven't fixed the warnings because another GSoC16 student is working on the package and because outside the scope of my project, but the problem can quickly be solved collaborating at the official end of this experience.