Planet Octave

May 13, 2017

Michele Ginesi


The expint function

The exponential integral

This is the first function I worked on in Octave: you can find the discussion here. Even if it will not be part of my GSoC project, I think it would be interesting to show how the function has been rewritten.

Definition and main properties

The canonical definition for the exponential integral is $$E_i (z) = -\int_{-z}^{+\infty} \frac{e^{-t}}{t}dt$$ but, for Matlab compatibility, expint compute the function $$E_1(z) =\int_{z}^{+\infty}\frac{e^{-t}}{t}dt.$$ The two definitions are related, for positive real values $z$, by $$E_1(-z) = -E_i(z) - \mathbb{i}\pi.$$ More in general, the function $E_n$ is defined as $$E_n(z) = \int_1^{+\infty}\frac{e^{-zt}}{t^n}$$ and the following recurrence relation holds true: $$E_{n+1}(z) = \frac{1}{n}\left(e^{-z}-zE_n(z)\right).$$ Moreover $$E_n(\bar{z}) = \overline{E_n(z)}.$$


To compute $E_1$, I implemented three different strategies:
  • Taylor series (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.11, p 229): $$E_1(z) = -\gamma -\log z -\sum_{n=1}^{\infty}\frac{(-1)^nz^n}{nn!}$$ where $\gamma\approx 0.57721$ is the Euler's constant.
  • Continued fractions (Abramowitz, Stegun, "Handbook of Mathematical Functions", formula 5.1.22, p 229): $$E_1(z) = e^{-z}\left(\frac{1}{z+}\frac{1}{1+}\frac{1}{z+}\frac{2}{1+}\frac{2}{z+}\cdots\right)$$ or, in a more explicit notation $$E_1(z)=e^{-z}\frac{1}{z+\frac{1}{1+\frac{1}{z+\frac{2}{1+\frac{2}{z+\ddots}}}}}$$ This formulation has been implemented using the modified Lentz algorithm ("Numerical recipes in Fortran 77" p.165).
  • Asymptotic series ("Selected Asymptotic Methods with Application to Magnetics and Antennas" formula A.10 p 161): $$E_1(z)\approx \frac{e^{-z}}{z}\sum_{n=0}^{N}\frac{(-1)^n n!}{z^n}.$$ A difficulty about this approximation is that the series is divergent, and that is the reason why the sum goes up to a certain $N$ "big" and not up to infinity.
Then I tested them in the square $[-100,100]\times[-100,100]$ of the complex plane, comparing the result with the ones given by the Octave symbolic package. With these informations, I identified three zones of the complex plane: in each zone one strategy is better than the others.
Then the final function divides the input into three parts, accordingly to the position of the same in the complex plane, and compute the function using the opportune strategy.

by Michele Ginesi ( at May 13, 2017 07:16 AM

May 10, 2017

Michele Ginesi


Introduction to the project


About me

My name is Michele Ginesi, I obtained a Bachelor's degree in Applied Mathematics in Verona, Italy. Now I am getting a Master degree in Mathematics in the same university.
I was selected to partecipate to the Google Summer of Code under GNU Octave for the project Make specfuns special again.

About the project

Special functions are an interesting and important topic in Mathematics, and so it is fundamental to have a way to compute them in an accurate way.
Some examples of special functions (whith some important application of the same) are:
  • Gamma function $\Gamma$. This is one of the most important: it can be viewed as the extension of the factorial function ($\Gamma(n)=(n-1)!$, for $n\in\mathbb{N}$) and it is a component in various distribution functions in probability theory.
  • Beta function $B$. This was the first known Scattering amplitude in String theory.
  • Bessel functions. These are the canonical solutions $u(r)$ of the Bessel's differential equation $$r^2\frac{d^2u}{dr^2}+r\frac{du}{dr}+(r^2-\alpha^2)u=0 $$ for any complex number $\alpha$. At a first approach may seems that they are an end in themselves, but actually the Bessel equation describes the radial component of the two dimensional wave equation ($u_{tt}-c^2\Delta u = 0$) applied on a disc.
The most common strategies used to approximate special functions are Taylor series and continued fractions, sometimes asymptotic series. Some particular functions can be implemented via recurrence formula or other type relations with other functions (for example, $B(z,w) = \Gamma(z)\Gamma(w)/\Gamma(z+w)$).
This project will be divided into three main parts:
  1. Fix the already known bugs related to special functions (e.g. #48316, #47800, #48036).
  2. When the known bugs will be fixed, I will proceed to add new tests to make sure that all the functions are accurate.
  3. Fix the new problems/bugs that, eventually, will be found during the second phase.
The main reference for this work will be Handbook of Mathematical Functions by Irene Stegun and Milton Abramowitz which contains all the functions I will work on completed with (almost) all the expansions that are needed to implement them; and Handbook of Continued Fractions for Special Functions by Annie Cuyt, Vigdis Brevik Petersen, Brigitte Verdonk, Haakon Waadeland and William B. Jones.
To test the functions I will use, in addition to the Handbook of Mathematical Functions, Sage and the Octave symbolic package to get a reference value.
Here you can find my repository on which I will work.

by Michele Ginesi ( at May 10, 2017 05:24 AM

May 01, 2017

Carnë Draug

My Free software activities in April 2017


Still reviewing the pending bugs and patches for the Image package:

  • fixed graythresh Otsu method when images have all the same value (bug #45333) by Avinoam Kalma.
  • new function imsharpen by Avinoam Kalma.
  • fixed normxcorr2 returning 0 for regions where image has all the same value (bug #50122) by Hartmut Gimpel.
  • fixed computation of ellipse properties (bug #49613) by Avinoam Kalma.
  • imreconstruct should clip marker values higher than mask instead of throwing an error (bug #48794)

Somehow got myself responsible for more things. Now also have access to Octave's server where the repos, wiki, and planet is hosted, to help out when needed.


Start packaging Algorithm::SVM (libalgorithm-svm-perl) which debian-med required to package psrtb.

Changed libbio-eutilities-perl package to remove unnecessary dependencies so it can be easily ported to older versions of Debian.


Fixed bug #781538 for correct highlighting of include type of LaTeX commands.

May 01, 2017 12:00 AM

April 01, 2017

Carnë Draug

My Free software activities in March 2017

This month I finished my PhD thesis


Attended OctConf 2017 at CERN where I presented how I use Octave at work for microscope image analysis. The slides are online although to be honest, I prefer to keep the slides with only a few images and talk through them wile waving my arms. I am unsure how useful they are outside the context of me talking.

Also at OctConf 2017, I gave two workshops: one on preparing Octave packages, and another on Autotools. The workshop on preparing Octave packages was a repeat of the same I gave at OctConf 2015. The workshop on Autotools was aimed to new Octave hackers and focused on how Octave makes use of Autotools.

With the thesis writing out of the way, I have started to review the pile of bug reports and patches to the Image package, submitted by Hartmut Gimpel and Avinoam Kalma over the last year:

  • fixed Image package to build against development Octave (bug #50180)
  • fixed bwperim handling of dimensions of length 1 (bug #50153)

April 01, 2017 12:00 AM

March 14, 2017

John W. Eaton

Timeline of Octave Releases

Someone asked me in private email why randi isn’t available in Octave. They were using 3.2.4 on Ubuntu 12.something and quoted a 6-year old email from the mailing list that said it was added in Octave 3.4. Maybe they are a time traveler from the past? Anyway, after fumbling around for the n-th time trying to find out when 3.2.4 was released (the answer is, a long time ago) I created this page to make those kinds of searches easier.

by jwe at March 14, 2017 07:25 PM

Silent Running, or, Back to the X-Term of the Future

I bought one of these:

put 8GB of memory and a 120GB SSD in it and mounted it on the wall.  So now I have (more or less) a modern-day quad-core X-Terminal with a gigabit network connection running Debian and Gnome.

Next, I moved my former noisy, heat-generating deskside system to the basement and, voila, sweet silence!  All I hear now is the residual ringing in my ears that is probably the result of attending too many loud rock-n-roll shows in my youth.  But that’s a story for a different time…

My home directory remains on the system that used to sit beside the desk.  It has a 6-core AMD CPU and a pair of disks in a RAID-1 configuration.  Accessing the files over the gigabit network seems fine, at least for normal usage.  And even though the new ZBOX has pretty modest specs, it is more than sufficient for most of the things I do: email, web browsing, and running Octave for testing.  I don’t use it  for building Octave.  For that I have a couple of other systems also located down in the basement where I can’t hear them.

Total cost was about $250.  I figured that wasn’t too much to pay for silence.  And if it turned out that the performance wasn’t great, I could always save up a few more pennies to buy a different model that has a faster processor and/or more memory and then reuse this box as a firewall – it does have 2 separate ethernet ports, as if it were intended to be used for that purpose.  Or maybe as an HTPC since it also has an HDMI output and Intel graphics that are supposed to support up to 4k graphics (I don’t know, my monitor works great with this system but it is just 1980 x 1200).

Setup was pretty easy and Debian installed without any problems from a USB stick.

It takes about 35 seconds to boot up to a gdm3 login screen.  That includes a few seconds each for POST and grub delays.  I could eliminate those and save a few seconds but then it would be harder to enter the BIOS or change boot parameters.  And how often do I reboot anyway?

From the login screen to a terminal window is another 15 seconds or so.  And that’s all with a 1.6GHz Celeron processor.  If I’d splurged and spent a few hundred more for the Core i5 version, I guess I’d be up and running in about 30 seconds?

by jwe at March 14, 2017 03:16 PM

March 01, 2017

Carnë Draug

My Free software activities in February 2017


New release of Pod::Weaver::Section::GenerateSection to make it work in Windows and newer perl versions.

New release of Pod::Weaver::Section::Legal::Complicated with a user testsuite. Until now it was an author only test because I'm afraid of changes in Dist::Zilla and Pod::Weaver but I think it's more important to have the tests available to cpantesters.

New release of Dist::Zilla::PluginBundle::BioPerl which drops deprecated plugins, adds a minimal testsuite, fixes an indexing problem in CPAN, and broken links for the bugtrackers.

Helped Pod::POM::View::Restructured use Dist::Zilla to handle their licence issues by making use of own [Legal::Complicated] plugin.


As part of Debian Perl group, I have packaged:

- libpod-weaver-section-legal-complicated-perl (I am upstream)
- libpod-weaver-section-generatesection-perl (I am upstream)
- libdist-zilla-role-pluginbundle-pluginremover-perl
- libconfig-mvp-slicer-perl
- libdist-zilla-config-slicer-perl

Took over libbio-eutilities-perl from the debian-med team since I am also upstream. They packaged it when I requested at a time I didn't know.

March 01, 2017 12:00 AM

February 20, 2017

John W. Eaton

August 30, 2016

Cristiano Dorigo

CGS, TFQMR and final revision of all the codes

Dear all,

During these weeks I worked on improving CGS, which I rewrite from scratch. Then I implemented TFQMR, which is not available in Octave. Finally I made a revision on the codes that I improved the past weeks.
All the file revisioned (gmres, pcg, bicg, bicgstab) and the new ones (cgs and tfqmr) are available here:

in the folder definitive_codes


As written in my second post, the actual Octave's CGS implementation is not very clear. Indeed I didn't understand most of the code (I was not able to recognize how the is was adapted). Moreover I tested it and the results are not correct (in my second post there is an example of this).
Then I decided to rewrite if from scratch following closely the algorithm given by Saad in the reference book suggested.
I wrote it following the pattern that I used in the other codes (i.e. to wrote three subscripts, each for the different cases: all matrices, all functions handle or mixed). Iit has also the same characteristics of the other methods improved, for example: the approximation returned is the one with the minimum residual, there are five types of flag (0 convergence, 1 maximum number of iterations reached, 2 preconditioner singular, 3 stagnation, 4 division by zero during the process)... I don't go in deep about this details to make less boring this post, since they are already explained in my previous posts. 


In my last post I wrote that the mentors told me to implement a new algorithm not available in Octave. Since I had no suggestion about which method write, and the community had no particular needs, I chose the transpose-free qmr algorithm because it is the only methods treated in the Saad book.
Also with this algorithm, I use the same pattern of the others (i.e. the three subscripts) and synchronized in such a way it has the same type of input and output as the others.

One mention to the number of iterations effectively performed. In the Saad book, the method has different behaviours if the iteration number is even or if it is odd, but at every iteration the approximation is computed. Making some test in Matlab with its tfmqr, I noticed that if the method converges, the vector of the residuals in long two times the number of iteration that Matlab tells that are necessary to convergence. Then I think that Matlab counts as one iteration the block (odd iteration + even iteration), i.e. I think that Matlab_it = odd_Saad_it + even_Saad_it. In this way the iterations performed by Matlab is two times the effective iterations performed. I decided to count the iterations in this way to be compatible.


Finally I checked (one more time) the codes that I improved in these weeks (i.e. pcg, gmres, bicg, bicgstab).
Indeed I updated the documentation (for example: making all similar, giving more details about the preconditioning, adding the reference book,...).
Then I checked the codes and I tried to make it more readable and to adapt it at the Octave conventions.
Then I added some tests in all the codes, mostly to check if all the subscripts works (i.e. checked every combination of A matrix/function, M1 matrix/function, M2 matrix/function), but also to check that the flags works well (I checked only the flag 0, 1 and 2, since the others are difficult to reproduce) and also to check that these methods solves also complex linear systems.

The mentors suggest me to make a patch with the methods improved (pcg, gmres, bicg, bicgstab and cgs), in such a way it can be easily the review of them and if it is in time (maybe) they can be included in Octave 4.2.
This patch is available here

and it is the file SoCiS16_Improve_iterative_methods.diff

I made this patch from the Octave changeset 22406:dc4e2203cd15 of the Octave 4.1.0+ version.


Since the official deadline for the SoCiS  is tomorrow (the 31-th of August), this is my last post in this blog.

To make a small summary, during this summer of codes I make an improving on
most of the iterative methods to solve linear systems. In particular I improved pcg, gmres, bicg, bicgstab and cgs.
I wrote for all of these methods three subscripts to make them more efficient, dividing the particular cases "matrix", "function handle" and "mixed", in such a way that, according to the input data, the method uses the appropriate subscripts to make computations.
I tried to make these methods more compatible with Matlab (for example noticing that Matlab gives as output the approximation with the minimum residual and not the last performed, as used in Octave, or adding the possibility of "flag =2", i.e. singular preconditioner).
I tried to fix some "non-clear" situation that I noticed, for example pcg does not noticed if the matrix was not Hermitian positive definite if it was a complex matrix, or to face the problem of a division by zero in bicg, bicgstab and cgs.
Then I implemented tfqmr which is available in Matlab, but not in Octave.

Unfortunately, there are other two methods that needs an improve (pcr and qmr), but due to unexpected problems with the above methods I had not enough time to improve also these. When my scholastic commitments are not so big, I want to try to improve also these two methods.

I want to thank all the Octave community that gives me this fantastic opportunity to face a big project as this and then to learn  much theoretically and practical things about these iterative solvers for linear systems.
I want to thank specially the mentors, who were always available to give suggestions, advices and ideas of how to proceed in the work, mostly when came some unexpected problems.
I hope that this improving will be useful to Octave, and to contact me if there is something not clear or wrong.
I also hope to continue to help the Octave community in the following times.

One more thanks to all of you.

by Cristiano Dorigo ( at August 30, 2016 10:34 AM

August 23, 2016

Francesco Faccio

Summary of the work - GSoC

MathJax TeX Test PageHello!

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:

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


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. 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.


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.

by Francesco Faccio ( at August 23, 2016 09:04 AM

Chiara Segala

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 ( 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/ : 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.

by Chiara Segala ( at August 23, 2016 03:09 AM

August 22, 2016

Amr Mohamed


GNU Octave – Geometry Package – GSoC 2016

This Repository ( ) is a clone(fork) of the Geometry Package which is a part of the Octave Forge Project (

This fork adds new functions to the official Geometry Package as part of GSoC (Google Summer of Code) 2016.

The official Geometry Package can be found here:

Added Functions

  • polyjoin: Convert cell arrays of polygons into NaN delimited coulmn vectors.
  • polysplit: Convert NaN separated polygon vectors to cell arrays of polygons.
  • ispolycw: Check whether the polygon(s) has(have) clockwise orientation or counterclockwise orientation.
  • poly2cw: Convert Polygons to clockwise contours(polygons).
  • poly2ccw: Convert Polygons to counterclockwise contours(polygons).
  • polybool: Perform Boolean Operations on polygons.
    Valid Operations are: (intersection,and,&), (union,or,|,+,plus), (subtraction,minus,-), (exclusive-or,xor).
  • poly2fv: Convert Polygons into Faces and Vertices Representation.


Community Bonding:

As mentioned in my proposal , I started working on the project 3 weeks before the official start of the GSoC program to compensate the three weeks that overlap with my final examinations.
I started my implementing the polysplit & polyjoin functions as they don’t depend on other functions.
These functions change the representation of the polygons from one form to another.
Since the polysplit & polyjoin functions don’t depend on Boost Geometry, they are implemented as .m files and placed in the inst/polygons2d/ directory.

First Coding Period:

During the first two weeks of this coding period, I continued the implementation of the polysplit/polyjoin functions adding help messages and regression tests.
Then, i worked on the first function that is implemented as a cc file which is the ispolycw function.

As all the package functions were implemented as .m files, i created a new directory src/ so that i can add .cc functions to the package.

The ispolycw function depends on calcualting the signed area of the polygon using the Boost Geometry function area(polygon).
The orientation of Self-intersecting polygons is defined as the respective orientation of the leftmost top point.If these points are co-linear ,  the signed area of the polygon is used to determine its orientation.

To compile the ispolycw function and generate .oct file that can be run by Octave, A new Makefile was added to the src/ directory.
The Makefile initially used the socket’s package implementation of the Makefile.

Two weeks after the start of the GSoC and we were already ahead of schedule. So we decided to start working on the polybool function as it is the main function of the whole project.
We decided to start by implementing the polybool function to support doing  boolean operation between only two polygons.

Second Coding Period:

By the start of the second coding period, we changed the implementation of the polybool  function to use the Boost Geometry multipolygon concept so that it can support applying boolean operations on multiple polygons with inner holes.

During the first three weeks the  second coding period, simple benchmarks were used to test the polybool  function.
And consequently multiple implementation bugs were fixed. One of these bugs was that a single inner ring was initially matched to multiple outer rings that surrounds this inner ring.
So a boolean array is currently used to ensure that each inner ring is matched to one and only one outer ring.

Then, i wrote two simple functions poly2cw poly2ccw which are used to toggle the orientation of polygons and thus convert inner rings (holes) into outer ones and vice versa.

Then in the fourth week of the second coding period, i decided to make better understanding of Makefiles , how to use them to build only changed files, how to organize dependencies between targets and how to use variables.
Also, i tried to use the dissolve extension of the Boost Geometry to handle self intersecting polygons.
It took me a few days to manage to use it as the Boost Geometry version on my system (Ubuntu 14.04) was 1.54 while the dissolve extension is developed based on Boost Geometry libray of version 1.60
After that, i adapted the ClipperLib benchmarks to Octave so that we can compare my implementation of these functions to the ClipperLib / other programs performance.

By the end of July, i  started working on the poly2fv function which is used to triangulate polygons so that they can be better visualized using the patch command.
This function returns a 2D matrix consisting of the triangles’ vertices and another 2D matrix representing the indexes of each triangle respectively.

For the first two weeks of August, i reviewed all the newly implemented package functions, added missing regression test , improved the help messages , added more input format checks , added warning messages when some boolean operations fail to be done mainly because of self-intersections that can’t be resolved.

Then during the last two weeks of the program, i had to learn how to use autotools so that i can use a configure script to test whether the Boost library is installed on the system, and check whether the dissolve function can be used to resolve self-intersecting polygons or not.
It was a great experience to manage to use the autotools to manage the generation of the final Makefile that is required to generate .oct files from the implemented .cc files.

Week Task Working hours / week
22-4 / 29-4 Understanding the organisation conventions – Prepare a complete plan for each function (Input/Output requirements – Optimal implementation based on memory/time metrics) 20 hour/week
30-4 / 6-5 Implementing the polysplit function and testing it. 20 hour/week
7-5 / 13-5 Implementing the polyjoin function and testing it. 14 hour/week
14-5 / 10-6 * Finalize the implementation of the conversion functions (polysplit-polyjoin) 7 hours /week
11-6 / 19-6 Implementing and testing the ispolycw function. * 40 hours/week
20-6 / 27-6 Midterm Evaluation week 40 hours/week
25-6 / 1-7 Working on the poly2cw and ispolycw scripts 40 hours/week
2-7 / 8-7 Performing severe testing to all the implemented functions 30 hours/week
9-7 / 15-7 Working on the polybool script. 30 hours/week
16-7 / 22-7 40 hours/week
23-7 / 29-7 Implementing and testing the poly2fv function. 40 hours/week
30-7 / 5-8 A 2 week buffer for finalizing the scripts , debugging them and writing the documentation. 40 hours/week
6-8 / 12-8 40 hours/week
13-8 / 23-8 Tidying the code, writing tests, improving the  documentation and submitting the code sample. 40 hours/week

* Final Examinations – Schedule isn’t announced yet – will inform the mentor whenever it is announced.


All the newly added functions has the following set of conventions:

  1. Outer Rings of the polygon has a clockwise orientation, while Inner Rings has a counterclockwise orientation.
  2. Description of the polygons is represented by two arguments X and Y. Each argument can be represented using Cell Arrays or nan Delimited Vectors.
  3. The functions don’t require that the last point is the same as the first one.

Example: To represent a square as shown in figure:

— NaN Delimited Vectors:

x=[0 0 3 3 NaN 1 2 2 1];

y=[0 3 3 0 NaN 1 1 2 2];

— Cell Arrays:

x={[0 0 3 3],[1 2 2 1]};

y={[0 3 3 0],[1 1 2 2]};


polybool & poly2fv functions use the Boost Libraries.

The recommended version is 1.60

Getting Started

To test the updates to the Geometry package:

  • Make sure you have Boost Libraries installed.
  • Clone the package
hg clone
  • Run ./bootstrap in the root directory of the package.
  • Create a tar file
tar -zcf geometry.tar.gz octave-geometry/
  • Run Octave.
  • In Octave Command Prompt, write
pkg install geometry.tar.gz
  • Load the package
pkg load geometry


All the benchmarks where run on a Ubuntu 16.04 computer with Boost Library 1.60.

All the scripts can be found here:

As my blog doesn’t support tables, Kindly check the Tests/Notes which can be found as a single pdf file here:


-poly2fv depends on polytri library which isn’t developed to throw exceptions when it can’t triangulate polygons so sometimes poly2fv calls may cause Octave to crash.

by amrkeleg at August 22, 2016 04:07 PM

Abhinav Tripathi


During Google Summer of Code 2016, I worked on the Symbolic Project under the organisation Octave. The project was directed towards adding a new IPC mechanism into symbolic – Pytave. During the course of GSoC work was done in 2 repositories (octsympy on github & pytave on bitbucket).

Code Changes:

All of the changes done by me can be viewed at the following links:

Summary of changes:

As we were to use pytave in Symbolic we worked on it to make it more mature. Add more functionalities to provide a better interface for Symbolic to deal with. Also many new features for conversion to and from python objects were introduced. The most important addition would be @pyobject and its functionalities.
On symbolic side, we used the features of pytave as much to our benefits as possible. We currently only have extra conversion for @sym objects which are specific to Symbolic. For rest of the octave types we use pytave‘s internal conversion.

Detailed PR wise changes:

To know about changes done before August, please read my other posts on this blog in the sequence: post-1, post-2, post-3 and post-4. Here I will summarise PR wise work done in August only.

In pytave, the PRs #33, #31 and #32 were merged which were under review in july. The PR 33 did some cleanup and added tests. PR 31 adapted existing tests since after inclusion of @pyobject we were now converting  python – lists, dicts and tuples; to @pyobject in octave. PR 32 fixed some doctests which were either not working in Python 3 or gave errors due to the latest @pyobject changes.
PR #24 added a new functionality to specify a namespace to be passed in py* functions to ensure that any changes in symbols/objects happens only in that namespace. This required some refactoring which was done in the same PR in a separate commit.
Then PR #40 added a method to allow comparison of 2 pyobjects. This would be useful in writing tests that need a desired output.

On the Symbolic side, PR #512 is under review which takes care of updating out existing methods to use the new pytave features. We have had a good amount of changes in pytave which now allow us to do most of the things without the need to execute and python code on our own (using pyexec). We now do it using py.* calls.

Goals not met:

Almost all the goals were met during the project. Only 1 major goal was not met which was building pytave on windows. Efforts were made towards that but later on it was more fruitfull to add new functionalities that would make the projects more mature. I would definetly try to work on this afterwards and hopefully in near future we will have fully working pytave on windows too.


by genuinelucifer at August 22, 2016 03:50 PM

Barbara Lócsi


GSoC 2016 is about to end so in this post I will sum up the work we had done during this period. I will write about 3 major topics: eig, gsvd and chol2inv. The original project was about eig, and when it was done we moved on to the other ones.


The goals of the project were:
Certain calling forms of the eig function are currently missing, including:
  •  preliminary balancing
  • computing left eigenvectors as a third output
  •  choosing among generalized eigenvalue algorithms
  •  choosing among return value formats of the eigenvalues (vector or matrix)

Calling forms:

We finished these goals before the end of GSoC, and most of the work is already merged to the main octave repository:
In my repository there are some other to be reviewed commits that reduces code duplication with the use of templates.
You can find some information about the individual tasks on this blog in previous posts, and also here:


The function gsvd was imported from Octave-Forge linear-algebra package. However it uses deprecated LAPACK functions: *ggsvd. Moreover it is not compatible with Matlab’s gsvd. For more information check these discussions:
  1. Deprecated LAPACK routines (GSVD)
  2. new gsvd function incompatible w/Matlab
  3. bug #48807: new gsvd function is incompatible w/Matlab
Work already merged to the main Octave repository:
Other commits that are not yet merged:


chol2inv couldn’t handle sparse matrixes (bug # 36437), we submitted a fix for this bug and tests for chol2inv.

GSoC 2016 with GNU Octave

This was an exeptional opportunity to gain experience and get into open source development and learn about Octave. I am really happy for this project and the time we spent working on it. I would like to thank for my mentors Carnë Draug and Nir Krakauer, without their help and guidance non of this would have been possible.

by Barbara Lócsi ( at August 22, 2016 02:02 PM

August 16, 2016

Francesco Faccio

How to replicate test


Here's a little guide to replicate the test I made for ode15i


Download last version of Sundials (only IDA module is required (IDA version > 2.8.0), but SundialsTB, idas and kinsol will be used for some test). Follow the install guide of Sundials to build IDA module with KLU (you will have to specify the lib and include directories where you have previously installed KLU). Use SUNDIALS_PRECISION as double (default) in configuration.


Clone my octave repo you can find here:
and build Octave from source. KLU is already a Octave dependence, but you can specify other versions you have installed at configuration time.


Here you can find the scripts and C functions I used to test ode15i.
You can test my oct-file running the scripts in folder Octave and the same scripts in Matlab folder for testing Matlab's ode15i.

To test the C examples of Sundials compile the C file you can find in folder C, linking against libklu, libsundials_ida and libsundials_nvecserial. If you installed these libraries in default path, type:

$ gcc -c idaHeat2D_klu.c idaRoberts_dns.c
$ gcc -o idaHeat2D_klu idaHeat2D_klu.o -lsundials_ida
-lsundials_nvecserial -lklu -lm
$ gcc -o idaRoberts_dns idaRoberts_dns.o -lsundials_ida
-lsundials_nvecserial -lklu -lm
$ ./idaRoberts_dns
$ ./idaHeat2D_klu

Sundials provides some mex interfaces for solving the benchmark problems. In order to test that, you have to edit file SundialsTB/install_STB.m changing mex compiler (near row 20):

mexcompiler = 'mkoctfile --mex';

Install Octave and run from Octave install_STB.m script. Choose to compile IDAS and kinsol interface and install the Toolbox. A startup file is generated where you chose to install SundialsTB: run it.

Now you could be able to run the script in MEX folder, which uses Sundials mex interface for solving the dense problem.

For problem ida_klu in Octave and Matlab folder, you can choose to test the solver with a sparse jacobian function, a dense jacobian function, a cell jacobian dense or a sparse jacobian dense.

Since I added the cell jacobian option in ode15i recently and I changed some parameters in the benchmark problems, I have some new results:

Problem Robertson eq (tspan from zero to 4e10, abstol = [1e-8, 1e-14,1e-6], reltol = 1e-4, dense methods):

Matlab: 0.41s
Octave: 0.06s
C impl: 0.002s
Mex impl: 0.12s

Without jacobian:

Matlab: 0.45s
Octave: 0.07s
Mex impl: 0.15s

Problem Heat2D (tspan from zero to 10.24, reltol = 1e-5, abstol = 1e-8, sparse methods):

With a sparse jacobian function:

Matlab: 0.34s
Octave: 0.13s
C impl: 0.005s

With a dense jacobian function:

Matlab: 0.44s
Octave: 0.15s

With a dense jacobian cell:

Matlab: 0.45s
Octave: 0.08s

With a sparse jacobian cell:

Matlab: 0.34s
Octave: 0.07s

Without a jacobian:

Matlab: 0.56s
Octave: 0.47s

It seems a quite goode result, considering that C and Mex implementations are specificly designed for solving these problems. Since I use a m-interface before calling my oct-file, there is a little overhead, but it will be negligible for bigger more complex problems.

For these test I configured Octave with -O2 flag and all the results refer to the code till this commit

by Francesco Faccio ( at August 16, 2016 01:43 PM

August 03, 2016

Amr Mohamed

Blog Post

Hi all,
I am so happy that our project is on the right track.
We have managed to implement the final function (poly2fv) on time.

This function perform Constrained Delaunay triangulation that i useful for plotting polygons with inner holes.
The function is implemented as a cc source file ,and it uses the polytri library which can be found here:

We have also adapted some of the ClipperLib Benchmarks to Octave and here is a screenshot for one of them:



by amrkeleg at August 03, 2016 04:16 PM

August 02, 2016

Cristiano Dorigo


Dear all,

In this last period I worked on the BICG and BICGSTAB algorithm.
You can find the codes in the folder bicg_bicgstab at

Since the main loop of both of them have more or less 30-40 lines of code, I wrote them from scratch, because in my opinion was easier instead to modify and to adapt the original code at the many situations that occurred.

The strategy that I used is the same of PCG and GMRES: for each algorithm I wrote three sub-scripts to subdivide the cases which: (a) the matrix and the preconditioners are all matrices, (b) all functions handle or (c) some are matrices and some functions handle.

As for PCG and GMRES, I tried to make these codes compatible with Matlab.

In my BICG and BICGSTAB codes there are not actually the documentation, but I'll add them as soon as possible

Main modifications on BICGSTAB and BICG

Many of the changes that I applied at the BICGSTAB and BICG were quite similar to those that I made in PCG and GMRES. I make a little list to summarize them: 

  • Now the approximation of the solution X returned is the one that has the minimum residual and not the last one computed by the algorithm. Also the output variable ITER contains the iteration in which X was computed, not the number of iterations performed.
  • In the original Octave BICGSTAB code the stagnation criterion is
    resvec(end) == resvec(end - 1).
    Instead in the BICG code it is: res0 <= res1 (where res1 is the actual residual and res0 the previous)
    I changed this criterion with:
    norm(x - x_pr) <= norm(x)*eps
    (where x_pr is the previous iteration approximation)
    This criterion is more efficient, since it checks the variation of the approximations and not the variation of the residual, and it seems more similar to the "stagnation" described in the Matlab help. Moreover in BICGSTAB (and in BICG) there are many fluctuations of the residuals (they are not decrescent in general), so it's better to check directly the approximations.
    Due to these fluctuactions, the BICG criterion is wrong, since in this manner the Octave BICG usually stops also if it is not necessary.
  • I added the "flag = 2" criterion, the one that Matlab describes as "preconditioner matrix ill-conditioned", that in Octave is skipped. The check that I implemented is for the singularity of the preconditioner matrix M, since seems that in Matlab this flag appears only in this situation.
Now I will focus on some situations regarding BICGSTAB:
  •  I added the "half iterations" as in Matlab: studying the Saad implementation, I noticed that the approximation X is made in two steps at every iteration, indeed it is defined as:
    x = x + alpha*p + omega*s
    Then to be compatible with Matlab (since p/alpha and s/omega are not related and computed in different situations), I set x = x + alpha*p (and its residual) as the half iteration, and x = x + alpha*p + omega*s as the complete iteration.
    In this way  the "real" total number of iteration performed (recalling that ITER contains the iteration which X was computed) are
    (length(RESVEC) - 1) / 2.   
  • There was a problem when the matrices A, M1, M2 and/or the right-hand-side b are complex, indeed many times the residual blow up and the method does not converge. In the Saad reference book (where it was taken the Octave BICGSTAB, I think), the method is formulated only for the real case. The only difference is that using this formulation, some inner product are inverted (e.g. it uses <v, w>, instead of <w, v>). In the real case there is no problem since the inner real product is symmetric, but in the complex not, indeed the result of <v, w> is the complex conjugate of <w, v>.
  • In the Matlab help I noticed that there is the possibility of "flag = 4". I tried to reproduce it, then I used the same input data in the Octave BICGSTAB: the approximation X was all composed by NaN, and also for the residuals. Then I study in deep the Octave code and the reference book. I discovered that ih this situation there was some diviosn by zero probably caused by the Non-Symmetric Lanzcos Biorthogonalization (that is used in both BICG and BICGSTAB).
    With this Lanczos algorithm is possible the situation of a "serious breakdown" (not going into deep details, but it happens when a certain inner product is zero, instead of a non-zero number), and there is not easy examples to reproduce it. The refernce book suggest two possibility to avoid this situation: to restart the Lanczos biortho, or to use a difficult variant called "Look-Ahead Lanczos".
    I tried to study the "Look-Ahead" but in the Saad book is only mentioned the general idea, and it is not simple to implement and to understand quickly. Then I tried the second strategy it is not efficient (an easy experiment to try this strategy is: run BICGSTAB to find the last iteration without NaN, then run another time the BICGSTAB with the output X as initial guess), indeed the method does not converge and the residual does not decrease.
    Then (since in Matlab is not mentioned if it uses or not the Look Ahead) in the BICGSTAB implementation that I wrote I don't use any strategy and I set "flag = 4" when the guilty inner product is zero.
    To see an example of this situation in the Octave BICGSTAB see the file example_bicgistab.m
Instead, during the implementation of BICG there was not (for now) strange situations. For both BICG and BICGSTAB I tried the same tests that are given in their Octave implementations and I add some test to try the "flag =2" and "flag = 4" situations.


I know that to fix PCG, GMRES, BICG and BICGSTAB I used a lot of time. This because in their implementation emerged a lot of strange situations that needed a deep revision of these algorithms and many times also a discussion (with the mentors or the community) of how to proceed to solve particular problems. Also the mentors told me that they expected less time/work to improve and make a revision of the codes already implemented.
They also told me to spend the remaning time of the SoCiS to implement some methods not present in Octave. For this, I want to ask you if there is one of the suggested algorithms (minres, symmlq, symmetr, tfqmr, bicgstabl or lsqr) that you prefer to implement.

As usual, if there is something wrong/not clear, please contact me via e-mail or write a comment.


by Cristiano Dorigo ( at August 02, 2016 06:52 AM

August 01, 2016

Abhinav Tripathi


Work in the past month has been distributed over symbolic and pytave side.

The work that was merged to symbolic in this time was mostly to convert objects from python to octave types. PR #474 and #491 worked in this regard.  PR #490 helped to make the tests compatible with Pytave IPC. Then we also added Pytave IPC to the build system (travis) with the PR #489.
Also, we worked to improve the passing of variables to python by using pytave. PR #465 has mostly achieved the goal. But it is failing on travis. The cause was finally detected in PR #502 to be the passing of struct to pytave. I couldn’t get much time to work on this as I was working on pytave. And hopefully when pytave matures then this error will be resolved in pytave itself.

In pytave, there were some bug fixes which caused octave to segfault. In the PR #12, the issue was addressed when an octave value could not be converted to python which threw an exception and crashed octave. In the PR #15, the issue where the exception text from Python could not be properly extracted was handled. This was later improved with the PR #30 for the cases of crash with Python 3. In PR #29, the issue was handled when a large number (integer) was being converted from python to octave then boost threw positive_overflow exception.

Also, there were some feature improvements too. PR #13, added the feature that checks for a function called via pycall into main namespace and builtins namespace too for existence. After pyobject was merged into pytave, most of the work was related to pyobject. PR #17 completed the functionality of returning pyobject from py* functions by adding the same mechanism to pycall which was implemented in pyeval. PR #16 added the feature to check length of a pyobject. PR #25 dropped the existing conversion mechanism which converted python list, tuples and dicts to octave matrices and structs. Now they return a pyobject instead.

There are a few PRs currently under review. PR #33 removes some of the unused code from pycall and adds a few BIST tests which showcase changes due to PR #25. PR #32 fixes some doctests which are also affected by the changes in the same PR. PR #31 aims at adapting all the methods of pyobject to the changes which dropped the conversion of lists, dicts and tuples. In PR #24, we are trying to add a new feature to allow users to specify the namespace in which pyexec/pyeval would run their code. Hence allowing isolation of different codes to have the same names/functions by running in different namespaces. Currently it is buggy and doesn’t actually work as expected. The symbols in a local namespace are also accessible from the global namespace. Moreover any 2 local namespace are sharing the same symbols. Currently it just works in the case where 1 global and 1 local namespaces are required. In such case, they are having separate symbol tables (dictionary actually as it is python).

During this time we had one setback. We tried a lot to build pytave on windows but failed. It has been stopped currently but we will resume it after GSoC ends. Tatsuro helped a lot in this. Thanks to him we have gotten really close to successfully building pytave on windows. The complete details about what was tried and which errors occurred can be found on my wiki page where I described Building Pytave on windows.

The project has been going smoothly because of a great participation from the mentors and octave community. Hopefully we will be able to have a working release of symbolic with pytave IPC (atleast for linux) before the end of GSoC.

by genuinelucifer at August 01, 2016 06:04 PM

Francesco Faccio

Performance test


As klu library is now linked to ode15i, I could test the performance of the code both in the dense and sparse case. I tested my oct-file against Matlab ode15i m-file, the C implementation of sundials examples and the m-file relying on the mex interface of sundials. Although the code still contains unefficient use of memory, I think that the results are pretty interesting.

Before going into details, I have to specify that C and mex implementations are written for solving these specific problems, so their code is highly optimized (specific Jacobian matrices are constructed directly inside functions required by IDA, while in my oct-file I receive a function that when evaluated it returns the Jacobian matrices I have to pass to IDA). In other words, solving a problem requires more time than solving the problem.

The first test regards the solution of Robertson chemical kinetics problem, in which the Jacobian is a 3x3 dense matrix. I found a solution for RelTol = 1e-4, AbsTol = [1e-8, 1e-14, 1e-6], tspan = [0, 0.4] and I got (in mean):

0.015s for the oct-file
0.13s for the mex-file
0.002s for the C-implementation
0.19s for Matlab ode15i.

The second test regards the solution of a 2D heat equation. In this problem the Jacobian was a 100x100 sparse matrix with 330 non zero elements. I tested with RelTol = 1e-5, AbsTol = 1e-8 and tspan =[0, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64, 1.28, 2.56, 5.12, 10.24] and got:

0.42s for the oct-file
0.005s for the C-implementation
0.45s for Matlab ode15i.

Unfortunately Sundials doesn't provide a mex interface for klu, so I couldn't test a sparse mex file version for this problem.

We can now make some considerations:
-Oct-file is 10 times faster than mex-file and this confirms that mex files cannot reach the performance of an equivalent (in this case, of a more general) otc-file
-Our solver is already faster than Matlab m-file
-The C implementation is the fastest, but code is designed to solve only a specific problem. We can try to get closer to this speed.

In order to improve the performance of the solver, the next step is to try to reduce looping and the use of memory and to write a DAE class, using functors and lambda functions to avoid the use of global pointers.

by Francesco Faccio ( at August 01, 2016 03:54 AM

July 26, 2016

Chiara Segala

odexprb34 solver

As I did with the odexprk23 solver from the Runge-Kutta family, I implemented the official exponential fourth order method with adaptive time stepping from Rosenbrock family described in my postAlso in this case I used the function expmv and I wrote a code for the stepper (exp_rosenbrock_34), in which there are three calls to expmv function. The first two for the calculation of the second and third stage
U_{n2} = u_n + \frac{1}{2} h_n \varphi_1( \frac{1}{2} h_n J_n) F(t_n,u_n) + \frac{1}{2} h_n \varphi_2(\frac{1}{2} h_n J_n) \frac{1}{2} h_n v_n ,
U_{n3} = u_n +  h_n \varphi_1( h_n J_n) [F(t_n,u_n) + D_{n2}] + h_n \varphi_2 ( h_n J_n) h_n v_n ,
[Atilde, eta] = (1/2*dt, Jn, [Fn, 1/2*dt*vn]);
X = myexpmv (1/2*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:, 2) = x(:) + X;

[Atilde, eta] = augmat (dt, Jn, [Fn+D(:, 2), dt*vn]);
X = myexpmv (dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:,3) = x(:) + X;

The third call to expmv is for both the solution to the next step, and the lower order solution for the estimation of the error

u_{n+1} = u_n + h_n \varphi_1(h_n J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n \dots
\dots + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] + h_n \varphi_4(h_n J_n) [ -48 D_{n2} + 12 D_{n3}] ,
\hat{u}_{n+1} =  u_n + h_n \varphi_1(h_n J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] ,
v = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3), -48*D(:, 2)+12*D(:, 3)];
w = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3)];
[Atilde, eta] = augmat (dt, Jn, v, w);
p = size (v, 2);
q = size (w, 2);
ep = [[zeros(p-1, 1); 1/eta; zeros(q, 1)], [zeros(p+q-1, 1); 1/eta]];
X = [eye(d, d), zeros(d, p+q)]*myexpmv(dt, Atilde, [sparse(d, 1), sparse(d, 1); ep]);
X1 = X(:, 1);
X2 = X(:, 2);
x_next = x(:) + X1;
x_est = x(:) + X2;

I wrote also a code for the solver (odexprb34), according to odepkg standards and  I tested the order of the method with the usual example described in my post.

by Chiara Segala ( at July 26, 2016 07:40 AM

July 18, 2016

Juan Pablo Carbajal

Vectores y matrices

En esta clase introducimos dos elementos fundamentales del calculo númerico, el vector y la matriz.

En la clase de Variables vimos que una variable puede guardar varios valores al mismo tiempo. Por ejemplo
u = [0,1,2,3,4,5]
Nota: los decimales en Octave se escriben con punto, no con coma. Es decir 1/2 se escribe 0.5 y no 0,5.

En el ejemplo, u es una variable que tiene 6 elementos o componentes. A estas variables, que organizan los valores de forma horizontal, las llamamos vectores fila. Las llamamos vectores fila porque los valores están organizados en una fila. Si ejecutamos el siguiente código:
size (u)

Octave nos dice:
ans =
1 6

OK, vamos por partes... Primero el comando size(u) nos devuelve el tamaño de u (size quiere decir tamaño en Inglés). Y la respuesta fue 1 6 ¿Como leemos la respuesta? 
El tamaño de un vector se indica de la forma "filas, columnas". En este caso 1 6 nos dice que u tiene 1 fila y 6 columnas.

Pregunta 1:
¿Cuantas filas y cuantas columnas tienen los siguientes vectores?
a = [0,1,2]
b = [0.1 0.5 -0.1]
c = [-3,2]
d = [83]
e = [0,1] 
tiene 1 fila y 2 columnas.

¿Cómo armamos un vector columna? Es decir, ¿Cómo organizamos valores de forma vertical?
v = [0;1;2;3]
Este es un vector columna, donde los valores se organizan verticalmente. Quiere decir que para crear un vector columna tenemos que separar las componentes con punto y coma.

Pregunta 2:
¿Cuantas filas y cuantas columnas tienen los siguientes vectores?
a = [0;1;2]
b = [0.1; 0.5; -0.1]
c = [-3;2]

Para acceder a la 2da componente del vector a = [0 -1 2] ejecutamos lo siguiente
y Octave responde
ans =

Supongamos que a(2) == -1 fue un error y deseamos cambiar ese valor por 1. Simplemente hacemos
a(2) = 1

Con toda esta información ya puedes armar tablas, o en términos matemáticos, una matriz. Una matriz no es otra cosa que la concatenación de filas o columnas ¿A alguien se le ocurre cómo armar una matriz? Por ejemplo ¿Cómo crearías la siguiente matriz en Octave?
0 1
1 0
0 1

Noten que uds. ya usan las matrices es su vida diaria, por ejemplo una tabla como esta:
 1. Literatura4
 2. Química 5
Puede entenderse como una matriz que tiene 2 columnas (Materia, Nota) y dos filas (cada materia).

¿Preguntas? Escriban en la sección discusión o en los comentarios más abajo.

En la próxima clase veremos como representar la posición de un objeto como un vector y su movimiento en el tiempo como una matriz.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM


¿Qué es una variable?

En lenguajes de programación, una variable es un símbolo (e.g. una letra:"x" o una palabra:"valor1" ) que se utiliza para guardar información.
En GNU Octave es muy fácil crear variables, simplemente
asignamos la información deseada a la variable. 

Por ejemplo:

x = 1

guarda el valor "1" en la variable llamada x. El signo = se utiliza para asignar valores, uno podría leer el código de arriba como "equis guarda el valor uno".

Pregunta 1: 
¿Cúal es el valor de la variable y en el siguiente código?

x = 1;
y = x + 1

Las variables en GNU Octave no solo pueden guardar un valor, sino que también pueden guardar listas de valores (vectores y matrices o tablas). 


x = [1,2,3]
y = "Hola Salta la linda!"

La variable x guarda un
vector fila
con componentes 1, 2 y 3. La variable y guarda un arreglo de caracteres, en otras palabras: un vector de letras! 
A estos objetos se los llama "strings" en Inglés y yo prefiero esta palabra que la palabra "arreglo", porque con esta última parece que estamos haciendo negocios sucios (chiste ;D).

Pregunta 2:
¿Cúal es la diferencia entre u y v en el siguiente código?

u = [0.1 -1 5];
v = [0.1; -1; 5];

Bueno, con estos introdujimos rápidamente variables en GNU Octave. Hagan sus preguntas! (y no se olviden de experimentar! También hay mucha ayuda en la web).

Hasta la próxima!

PS: No se olviden que todo el código de las clases se puede ejecutar en una sesión de GNU Octave.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

¡Animarse a animar!

Todos disfrutamos de una buena animación, aquí veremos como podemos animar nuestros datos en GNU Octave. Les presento algunos ejemplos a modo de motivación:

Para lograrlo vamos a conocer la función comet y, para hacer nuestras propias animaciones, vamos mirar un poco en detalle los objetos generados por la función plot.

La función comet

En su forma más intuitiva, esta función toma como entradas tres argumentos: 2 vectores (o listas de valores) y una valor escalar (un número). Los dos vectores son los pares de valores (x,y) que se irán dibujando uno a uno. El argumento escalar define el tiempo de espera antes de mostrar el siguiente punto. Probemos lo siguiente:

T   = 10;
fps = 25;
t   = linspace (0, T, fps*T).';

l = 0.5;
w = 2*pi*0.5;

R = exp (-l * t);
x = R .* cos (w * t);
y = R .* sin (w * t);

fig = figure (1);
set (fig, "Name", "Espiral");

comet (x, y, 1/fps);

Si ejecutas el código, deberías ver algo como lo que se muestran en el siguiente video:
Pregunta 1:
¿Entiendes todo el código?

Anatomía de una animación
Para realizar cualquier tipo de gráfico en Octave necesitamos una figura (la ventana donde vemos el gráfico, creada con la función figure). Esta figura contiene un par de ejes coordenados (un "hijo", que se guarda en el campo de la figura llamado, en inglés, "children"); que en Octave se llaman axes y se que crean con la función homónima. Dentro de estos ejes podemos tener objetos gráficos, como líneas, puntos y polígonos.

Como se darán cuenta, las figuras pueden tener mas que un hijo, es decir una figura puede contener varios ejes permitiendo poner varios gráficos en la misma figura.

Resumiendo: una figura contiene un par (o más) de ejes coordenados, los cuales a su vez contienen objetos gráficos
Veamos un ejemplo

close all
fig = figure (1, "Name", "Figura 1");  # Figura con nombre
ax  = axes ();                         # fig adopta ejes automáticamente.

x = randn (10,1);                      # Graficamos algo
h = plot (ax, x, '.');                 # y lo ponemos en los ejes creados

Por supuesto, todo esto (excepto el nombre de la figura) se hace automáticamente simplemente ejecutando plot (x,'.') ¡Verifícalo!
En la última línea del código, le pedimos a la función plot que nos devuelva el objeto gráfico. Podemos ver las propiedades de este objeto ejecutando get(h) y obtendremos una lista enorme de propiedades del objeto gráficos y sus valores actuales. Fíjate que en esa lista hay campos con los nombres ydata y xdata. Estos campos contienen los puntos que se muestran en el plot. Si modificamos estos campos, modificaremos el gráfico. Por ejemplo podemos desplazar los valores un lugar hacia la derecha haciendo

t = get (h, "xdata");
set (h, "xdata", shift (t, -1));

Pregunta 2
¿Puedes modificar el código para realizar la animación que se muestra abajo (click para ver animación)?
gif animado

Hemos visto los principios de animación: prepara el gráfico y luego modifica su contenido.
En general es mala idea rehacer el gráfico a cada instante, pues es muy lento y las animaciones quedan mal.
Puedes animar cualquier propiedad de los objetos gráficos, ejes y figuras.

Si tienes alguna idea y necesitas ayuda no dudes en preguntar en el foro

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

¿Para qué me sirven los vectores?

En la clase de vectores comentaba como se puede crear un vector en GNU Octave, pero quedó en el tintero la motivación para el uso de vectores ¿Para qué me pueden servir esas listas de numeritos que llamamos vectores?

Cartas de amor...

Supongo que alguna vez han mandado una carta por correo postal. Para que la carta llegue a destino necesitamos proveer cierta información además del nombre del destinatario. Supongamos que la carta la enviamos dentro de Argentina. Necesitamos indicar la provincia a donde va la carta, la cuidad, la calle y el número de casa. En principio podríamos organizar la información en de esta forma

carta = [ "Buenos Aires","Lomas de Zamora", "Ceferino Namuncura", "150" ];

Ok, esto ya se parece a un vector pero las componentes son strings en vez de números. No hay drama, simplemente podemos indicar la provincia con un número del 1 al 23 (Argentina tiene 23 provincias) y hacer algo parecido con las ciudades y las calles. El número de casa lo podemos usar directamente. Es decir que el destino de nuestra carta podría representarse con el vector

carta = [ 23, 1352, 12345, 150 ];

donde cada componente fue reemplazada por un número según una tabla especificada. La primera componente del vector nos indica la provincia, la segunda nos muestra las ciudad, la cuarta nos indica la calle y la quinta componente la altura de la calle ¡Estas cartas viven en un espacio de 5 dimensiones!
Un conjunto de estas cartas se puede organizar en una lista de vectores, una matriz:

cartas = [ 23, 1352, 12345, 150; ...
          5, 130, 4, 756; ...
          12, 7, 2341, 29 ];  

Cada vector fila de esa matriz representa una carta. La primera columna de la matriz nos muestra las provincias, la segunda columna nos muestra las ciudades, la cuarta nos indica la calle y la quinta columna la altura de la calle. 

Panadero, quiero pan!

Ok, quizás es ejemplo fue muy abstracto. El siguiente ejemplo me lo pasaron por Facebook, también es abstracto, pero quizás es más fácil de entender. 

Pensemos por un momento en recetas de cocina, en particular en la lista de ingredientes de cada receta. En la tabla puse dos ejemplos de ingredientes para hacer pan

 Pan de cazuela  Pan fuerte
 1.Harina (gr)
 250 300
 2.Sal (gr) 5 20
 3.Levadura (gr) 5 20
 4.Aceite (ml) 0 50
 5.Agua (ml) 175 50
 6.Miel (ml) 0 100
 7.Huevos 0 0

Esta tabla ya nos muestra una representación de estas recetas como vectores. En este caso son dos vectores columna que en GNU Octave podríamos anotar como

recetas = [ 250 300; 5 20; 5 20; 0 50; 175 50; 0 100; 0 0 ];

Los ingredientes del pan de cazuela son la primera columna de esta matriz y los del pan fuerte son la segunda. Noten que estos vectores son de 7 dimensiones (7 ingredientes), pero vale preguntarse si la dimensión Huevos tiene sentido en el espacio de recetas de pan... ¿Hay recetas de pan que utilizan huevos?

Los dos ejemplos anteriores son bastante abstractos y si los has entendido me puedo imaginar que ya tienes muchas ideas de como representar otras cosas como vectores o matrices. 

Trayectorias espaciales

El siguiente ejemplo es paradigmático en el campo del cálculo numérico, y es muy relevante para este curso: representar la posición de objectos en el espacio!

Imagina una hormiga caminando sobre una hoja de papel A4 (ancho: 21cm, largo: 29.7cm). Para nuestra conveniencia hemos marcado una de las esquinas de la hoja como punto de referencia. En cualquier instante podemos leer la distancia entre la hormiga y la esquina de la hoja en las dos dimensiones del papel: ancho y largo. Ponemos la hormiga en un punto del papel y una vez por segundo miramos donde se encuentra. Podríamos observar algo así:

trayectoria hormiga

Figura 1: Posiciones de la hormiga en el papel. Cada punto representa la posición de la hormiga en las coordenadas del papel observadas cada 1 segundo. La línea punteada es la trayectoria de la hormiga observada con más frecuencia.

La posición luego del primer segundo la guardamos en una matriz con una fila y con dos columnas (un vector fila!)

p(1,1:2) = [ 10.7 14.8 ];

Luego de otro segundo obtenemos la nueva posición y la guardamos en la segunda fila de la matriz:

p(2,1:2) = [ 11.1 14.8 ];

Y así sucesivamente. Lego de 10 segundos, la matriz tiene 10 filas y 2 columnas. Esta matriz representa la trayectoria de la hormiga a intervalos de 1 segundo. Cada vector fila nos indica la posición en un dado instante de tiempo.

p = [ 10.7 14.8; ...
      11.1 14.8; ...
      11.4 15.1; ...
      11.2 15.0; ...
      11.0 15.3; ...
      11.8 15.8; ...
      11.6 16.8; ...
      11.8 17.0; ...
      11.0 17.3; ...
      10.8 17.2];

Creo que estos tres ejemplos deberían entusiasmarte para pensar qué cosas pueden representarse utilizando vectores y matrices y cómo hacerlo. Vale la pena preguntarse como podemos utilizar estas representaciones. Por ejemplo, usando la matriz de la hormiga ¿Cuál es la distancia entre la primera y la ultima posición de la hormiga?

¿Se te ocurren otros ejemplos? ¡Anótalos en los comentarios!

Nota: El archivo adjunto es un script de GNU Octave que genera la trayectoria de la hormiga y crea la figura que se muestra en esta clase.

by Juan Pablo Carbajal ( at July 18, 2016 07:57 PM

Funciones para generar valores

GNU Octave provee un gran número de funciones que son útiles para resolver nuestros problemas numéricos.
No tiene sentido enumerar aquí al enormidad de funciones de Octave, para eso está el manual y la ayuda interactiva que describimos en una clase pasada.  Los invito a que usen el foro de discusión para hacer preguntas del tipo: ¿Con qué función puedo hacer XXXX?

En nuestras clases irán apareciendo funciones, cuya funcionalidad la explicaremos en su debido momento (en realidad una vez que saben el nombre de la función pueden ir al manual o utilizar la ayuda interactiva).
En esta clase, introduciré dos funciones que nos permiten generar valores que incrementan de forma lineal o logarítmica.

Valores igualmente espaciados

La función linspace toma tres argumentos de entrada:

x = linspace (inicio, final, numero);

El argumento inicio indica el primer valor y el argumento final el último. El argumento numero indica cuantos valores generar entre inicio y final, todos ellos igualmente espaciados.

valores igualmente espaciados

Pregunta 1:
Si quiero valores que van desde 0 a 16 con un espacio entre ellos de 2, puedo ejecutar el código:

x = 0:2:16;

¿Cómo hay que llamar a linspace para generar los mismos valores?  ¿Puedes explicar la solución?

Valores logarítmicos

La función logspace funciona igual que la función linspace pero en vez de distribuir valores igualmente espaciados, distribuye los exponentes de los valores. Por ejemplo

x = logpspace (-1, 3, 5);

genera los valores 

x =

   1.0e-01   1.0e+00   1.0e+01   1.0e+02   1.0e+03

Si observan bien, pueden notar que es el exponente el que está igualmente espaciado (-1,0,1,2,3). Es decir, x(1) es 1x10-1 mientras que x(end) es 1x103.

valores espaciados logaritmicamente
Pregunta 2:
¿Qué valores obtengo si ejecuto el siguiente código?

x = log (logspace(-3,3,10))

Pregunta 3:
¿Cúal es el comando que genera valores logarítmicos entre 0.5 y 125?

No duden en dejar sus preguntas, sugerencias y comentarios en el foro!

¡Hasta la próxima!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¡Ayuda! ¡Socorro!

GNU Octave tiene un sistema de ayuda integrado. En esta clase veremos como usar el sistema de ayuda y otras fuentes de información. La mayoría de esta ayuda está en Inglés, y es importante saber leer este idioma a la hora de programar. Sin embargo muchas de las personas que trabajan en GNU Octave hablan español y están dispuestos a ayudar!

En la terminal

Mientras trabajamos en Octave hay varias formas de pedir ayuda sobre una función. Cada forma se aplica en diferentes casos.

Si sabes el nombre de la función tu mejor amigo es el comando help. Supongamos que sabemos que la función que queremos usar es la arco tangente (la inversa de la función tangente). En Octave, hay dos funciones que calculan la arco tangente: atan y atan2. Veamos la ayuda de cada función (en Octave 3.8.2):

> help atan
'atan' is a built-in function from the file libinterp/corefcn/

 -- Mapping Function: atan (X)
     Compute the inverse tangent in radians for each element of X.

     See also: tan, atand.

> help atan2
'atan2' is a built-in function from the file libinterp/corefcn/

 -- Mapping Function: atan2 (Y, X)
     Compute atan (Y / X) for corresponding elements of Y and X.  Signal
     an error if Y and X do not match in size and orientation.

     See also: tan, tand, tanh, atanh.

¿Entiendes la diferencia entre estas dos funciones?

A medida que vayas acumulando más experiencia con Octave, podrás recordar los nombres de la funciones que te resultan más útiles. Sin embargo puede ser que a veces no te acuerdes el nombre completo de la función pero solo el principio. En este caso podemos utilizar el autocompletar para obtener el nombre de las funciones. Supongamos que nos acordamos que una función empezaba con "sum" pero no recordamos el resto. en la terminal de Octave ingresamos este texto y luego presionamos dos veces seguidas la tecla TAB. El resultado es el siguiente

> sum
sum     sumblk  summer  sumsq

Octave nos da una lista de todas las funciones que empiezan con "sum" ¡Esto es de gran ayuda para encontrar nombres de funciones!

Otra situación bastante común es cuando no sabemos el nombre de la función pero sabes que el texto de ayuda contiene una frase/palabra clave que si conocemos. En este caso una herramienta muy útil es lookfor (significa "buscar" en Inglés). Consideremos un ejemplo. Supongamos que queremos buscar si Octave tiene una función para calcular la suma de una lista de números. Como frase a buscar (en Inglés) elegimos "sum of" (suma de) y ejecutamos:

> lookfor "sum of"
trace    Compute the trace of A, the sum of the elements along the main diagonal.
hypot    Compute the element-by-element square root of the sum of the squares of X and Y.
cumsum   Cumulative sum of elements along dimension DIM.
sum      Sum of elements along dimension DIM.
sumsq    Sum of squares of elements along dimension DIM.
md5sum   Calculate the MD5 sum of the file FILE.

Octave nos responde con una lista de nombres de funciones y la primera linea de la ayuda correspondiente. Una vez que tenemos los nombres de las funciones podemos utilizar el comando help para investigar en mas detalle que es lo que hace cada función.

Finalmente, puedes acceder al manual de Octave en la terminal utilizando el comando doc.

En internet

Como siempre estar conectado hace la vida más fácil. En internet podemos utilizar el manual online de octave para buscar ayuda sobre las diferentes funciones y funcionalidades. Es un buen lugar para consultar y puedes descargar la versión en PDF si acaso no tienes conexión todo el tiempo. He puesto algunas referencias en castellano en la sección de lecturas recomendadas. También puedes visitar el Wiki de octave.

Otro recurso excepcional es el chat room #octave en Aquí puedes charlas con otros usuarios y desarrolladores de todas partes del mundo. No es raro ver conversaciones en castellano, italiano y alemán. Puedes accerder al chat room utilizando un sitio web o cualquier cliente de IRC (por ejemplo Pidgin).

Finalmente, si prefieres conversaciones por e-mail puede utilizar la lista de mails

¡Utilizando todos estos recursos no puedes sentirte solo a la hora de aprender GNU Octave!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Repitiendo instrucciones: lazos "for"

Seguimos con el tema solicitado por Matías: control de flujo. Ahora vamos a ver como repetir un grupo de instrucciones utilizado el lazo for (yo les digo "loops", no lazos).
La idea es sencilla: uno tiene que repetir una acción muchas veces y no quiere escribir el código correspondiente a cada una de ellas.
Por ejemplo, supongamos que queremos sumar todos los valores en un vector:

x = rand (100,1)

El vector x contiene 100 números entre 1 y 0 elegidos al azar (ver ayuda de la función rand y preguntar en el foro si algo no se entiende). El objetivo es sumar todos estos valores. La forma ingenua de hacer esto sería escribir las 100 líneas:

y = x(1);
y += x(2);
y += x(3);

Y así hasta sumar todo los números. Está claro que esto no nos lleva muy lejos. Los lazos nos permiten repetir un grupo de comandos sin tener que escribir todas estas líneas. Utilizando el loop for podemos hacer la suma:

x = rand (100,1);
y = 0;
for i = 1:100
    y += x(i);

El loop define una variable de iteración i e indica que esta variable va desde 1 a 100. Esta variable se utiliza luego para indexar el vector x y acceder a cada uno de sus valores. Para vectores de otros tamaños el programa es casi el mismo

x = rand (500,1);
y = 0;
for i =1:500
    y += x(i);

Pregunta 1:
¿Cómo podemos escribir el loop para que funcione con cualquier tamaño del vector x?

El precio de la iteración

GNU Octave se basa en un paradigma de programación vectorial. Esto quiere decir que operaciones entre vectores y matrices son veloces y las iteraciones explicitas, como el loop for, son lentas. La idea es que en vez de hacer una iteración, podemos hacer operaciones muy rápidas utilizando más memoria. Es decir, cambiamos  iteraciones por espacio en memoria.

El ejemplo anterior se puede escribir de manera vectorial de la siguiente manera:

x = rand (500,1);
y = sum (x);

Este ejemplo en particular no cuesta más memoria que el loop porque en ambos casos construimos el vector x. Sin embargo, el loop podría ahorrar memoria si lo escribiésemos así:

y = 0;
for i =1:500
    y += rand ();

En este caso nunca construimos el vector x y ahorramos esa memoria. Esto podría ser útil en el caso que necesitas una cantidad tan grande de números aleatorios que no entran en la memoria de tu PC (si tienes 1Gb de memoria puedes tener alrededor de 1 millón de estos números). Sin embargo esto tiene un precio en tiempo y lo podemos verificar utilizando las funciones tic y toc. El siguiente programa compara los tiempos de ejecución de las dos formas de hacer esta suma:

N = 1000;
disp ('Loop')
y = 0;
for i =1:1000
    y += rand ();
disp ('Vectorizado')
x = rand (1000,1);
y = sum (x);

En mi PC esto da 

Elapsed time is 0.0164838 seconds.
Elapsed time is 7.60555e-05 seconds.

¿Cuánto da en la tuya? Puedes descargar el script al pié de la página.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¡Por favor Sr. Kirchhof, resuelva mi circuito!

En esta entrega vamos a utilizar lo que discutimos sobre ecuaciones lineales para resolver problemas con circuitos eléctricos.
En la figura siguiente vemos un circuito compuesto de 4 resistores y una fuente de voltaje.
El dibujo de la izquierda muestra el circuito con la fuente conectada al nodo 1 y el nodo 4 conectado a tierra. El dibujo de la derecha muestra un circuito equivalente. Nosotros podemos controlar la fuente y por lo tanto conocemos el voltaje Vi entre los nodos 1 y 4. Todos las otras 4 caídas de voltaje V12, V13, ... son incógnitas que tenemos que averiguar. También son desconocidas las 5 corrientes: la corriente I entregada por la fuente y  la corriente circulando por los 4 resistores I12, I13, .... Planteemos las ecuaciones de este circuito y resolvamos con Octave.

Ley Kirchhof del voltaje

La ley de Kirchof del voltaje (que es una ley que se puede derivar de las ecuaciones del electromagnetismo), dice que en ausencia de campos magnéticos, la suma de las caídas de voltaje en un circuito cerrado (en un lazo, en un loop) debe ser 0.
En la figura de arriba, a la derecha, he marcado dos lazos con flechas rojas. La dirección de las flechas es arbitraria y simplemente indica en que dirección el voltaje cae. Lo único importante es que una vez que hemos dado un sentido a las flechas no lo podemos cambiar (si lo cambiamos mientras escribimos las ecuaciones, estas van a estar mal). Es decir que si V12 es positivo, quiere decir que el nodo 1 tiene un voltaje más alto que el nodo 2. Si este voltaje es negativo, entonces es al revés, el nodo 2 tiene un voltaje mas alto que el nodo 1.
En el dibujo de la izquierda, marque otro lazo que no involucra a la fuente.

Escribamos las ecuaciones:

Las dos primeras corresponden a los lazos en la figura de la derecha. El primer lazo empieza en el nodo 1, va al 2 luego al 4 y vuelve al 1.  Cuando circulo en la misma dirección que la flecha el voltaje aparece en la ecuación con signo positivo, y cuando circulo en contra del sentido de la flecha, el voltaje aparece con signo negativo (¡voy en contra!). En el caso del segundo lazo (la segunda ecuación) me paseo por los nodos 4,3,1 y de regreso  4. al tercera ecuación es dar una vuelta pasando por todos los resistores. 
¿Se dan cuenta que la suma de las dos primeras ecuaciones genera la tercera? Esto quiere decir que la tercera ecuación no nos brinda nueva información: sabiendo las dos primeras ecuaciones, ya sabemos lo que dice la tercera. Por lo tanto la tercera ecuación no la vamos a utilizar (podrían no utilizar cualquier otra).
Aplicando nuestra habilidad para transformar sistemas de ecuaciones lineales en productos de matrices podemos obtener:

En este sistema parecen las 4 incógnitas de voltaje. En Octave podemos escribir la matriz de factores de la ley de Kirchhof de voltajes (LKV), de la siguiente manera:

LKV = [1 0 1 0; 0 1 0 1];

Ya tenemos las leyes de Kirchhof de voltaje, ahora veamos las corrientes.

Ley Kirchhof de las corrientes

La idea detrás de esta ley es que la corriente se conserva. Es decir que cuando corrientes se unen o se separan en un nodo, la suma se tiene que conservar. Para escribir estas ecuaciones nos paramos en un nodo y vemos que corrientes llegan al nodo (positivas) y que corrientes se van del nodo (negativas). En este caso también obedecemos el sentido de las flechas que definimos al principio. Repito, es importante mantener el sentido de las flechas hasta el final, pero el sentido es arbitrario. Una vez que elijen un sentido para la flecha, no lo cambien hasta que terminen.
Hagamos un ejemplo: la corriente que llega al nodo 1 viene directo de la fuente (corriente I); la corriente I13 también llega al nodo, mientras que la corriente I12 se va del nodo. Todo esto nos da la primera ecuación del siguiente sistema:

Las demás ecuaciones las pueden escribir siguiendo el mismo método para cada nodo. La ultima ecuación, nodo 4, también se puede deducir de las tres primeras (¿Pueden explicar cómo?), por lo tanto no la vamos a utilizar. El sistema escrito en forma matricial es:

Noten que tenemos 3 ecuaciones y 5 incógnitas de corriente.

En Octave podemos cargar la matriz de factores para la ley de Kirchhof de las corrientes (LKC) con el código:

LKC = [-1 1 0 0 1; 1 0 -1 0 0; 0 -1 0 1 0];

Hasta aquí hemos coleccionado 5 ecuaciones: 2 de LKV y 3 de LKC. Sin embargo tenemos 9 incógnitas, 4 voltajes y 5 corrientes. Para poder encontrar tantas ecuaciones como incógnitas necesitamos recurrir a la relación entre voltaje y corriente determinada por los resistores.

Relación entre corriente y voltaje

En el k-ésimo resistor, la relación entre la caída de voltaje en el resistor y la corriente fluyendo a través de el, viene dada por:

Como ven es también una relación lineal: la caída de voltaje es proporcional a la corriente. Escribiendo esta ecuación para cada resistor, obtenemos el siguiente sistema de ecuaciones:

Para construir esta matriz en Octave, es necesario dar el valor de resistencia de los resistores. Aquí elijo cualquier valor, uds. pueden usar otros:

R = [1 1 3 3];
VIR = [eye(4,4) diag(-R) zeros(4,1)];

En la última linea de comando utilizo 3 funciones que no hemos discutido aún. Para entender que hacen tenemos que mirar un poco la forma de las matrices en el sistema de ecuaciones. Si se fijan en la primera parte de la matriz podrán ver que hay unos (1) formando una diagonal. En la segunda parte también tenemos una diagonal formada por el negativo de las resistencias. La última columna de la matriz está llena de ceros. Ok, a las funciones:

eye(n,n): Esta función genera una matriz con n filas y n columnas que tiene la diagonal llena con unos.

Pregunta 1:
¿Qué matriz se genera si llamamos a esta función con dos argumentos distintos? Por ejemplo eye(2,3).

diag(v): Esta función genera una matriz cuya diagonal se llena con los elementos del vector v.

Pregunta 2:
¿Podés generar la matriz eye(4) con la función diag()?

zeros(n,m): Esta función genera una matriz con n filas y m columnas llenas de ceros.

Ensamblando todo el sistema

Quizás ya se han dado cuenta que ahora tenemos tantas ecuaciones como incógnitas. Ahora podemos armar el sistema de ecuaciones completo

Con color resalté cada uno de los bloques de la matriz y su procedencia. Los bloques marcados con un cero grande, indican que allí solo hay ceros.
En Octave podemos ensamblar esta matriz y resolver el sistema (¡Ojoi! antes de resolver necesitan darle un valor al voltaje Vde la fuente) con el siguiente código:

A   = [LKV zeros(2,5); zeros(3,4) LKC; VIR];
neq = size (A,1); # Numero de ecuaciones
b   = [-Vi; Vi; zeros(neq-2,1)];

# Resolver el sistema
x  = A \ b;
V  = x(1:4);
I  = x(5:8);
Ii = x(9);

disp ("Voltages"); disp (V); disp ("Corrientes"); disp (I); disp ("Corriente de la fuente"); disp (Ii)

Noten como he vuelto a utilizar la función zeros() para rellenar los bloques con ceros. La última linea de este código muestra los resultados en pantalla ¿Coinciden con tus cálculos manuales?

El archivo adjunto contiene todo el código aquí utilizado.

¿Encontraste un error en este post? ¿Tienes preguntas? Utiliza el foro de discusión.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Ecuaciones lineales? ¡No problem!

Uno de los problemas más recurrentes en el cálculo numérico es el de resolver sistemas de ecuaciones lineales. ¿Qué es una ecuación lineal? Si no conocen este nombre todavía, no se preocupen. Vamos a ilustrar todo desde el principio y de manera básica. También veremos una aplicación fundamental de las matrices.

La ecuación lineal

Los problemas lineales son muy comunes en la vida cotidiana. Son tan comunes que en el colegio nos enseñan a resolverlos con una regla: la famosa, y a veces odiada, regla de tres simple.
Esta regla resume la facultad de resolver proporciones. Ejemplo:

Problema 1:
Si un paquete de harina pesa 10 kg cuando está lleno ¿Cuánto pesa cuando solo queda un tercio del paquete?

En cada generación y en cada escuela se enseñan distintas formas de plantear este problema, pero la idea fundamental es que el peso del paquete es proporcional a la cantidad de harina. Esto quiere decir que mientras más(menos) harina, más(menos) pesa el paquete. Pero no solo implica que las dos magnitudes aumentan juntas, sino también que aumentan juntas de una manera particular: cualquier aumento de harina implica un aumento de igual proporción en el peso del paquete. Esto es lo que se llama una relación lineal

Donde P es el peso del paquete y V es la fracción de llenado del paquete (1=lleno, 0=vacío). El factor k es la proporcionalidad entre P y V. En nuestro problema la proporción es:

Otros ejemplos de proporcionalidades son la relación entre el volumen y la masa de un líquido (la proporcionalidad se llama densidad del líquido); entre el volumen y la altura de un cilindro (la proporcionalidad está dada por el área de la base del cilindro); entre el la altura de un fajo de billetes y la cantidad de billetes (la proporcionalidad es el espesor de cada billete), etc... seguro pueden imaginarse en un montón de relaciones lineales.

Una relación lineal que aparece en física del colegio secundario es la relación entre distancia recorrida en linea recta y rapidez.

Donde x es la posición final, x0 la posición inicial, v es la rapidez y Δt el tiempo transcurrido. Como ven, la distancia recorrida es proporcional al tiempo transcurrido; también podemos decir que la distancia recorrida es proporcional a la rapidez.

Sistemas de ecuaciones lineales

Consideren el siguiente problema.

Problema 2:
Encuentra los números tales que: si sumamos el doble del primero con la mitad del segundo, obtenemos 10 y si restamos el primero al segundo obtenemos 5.

Escribamos este problema en ecuaciones

Los números que estamos buscando son x,y  que al momento representan incógnitas. Podemos organizar las ecuaciones de la siguiente manera:

donde el punto indica producto escalar entre el vector fila y el vector columna (ver clase sobre el producto escalar). Ok, si miran estas ecuaciones con fuerza, se darán cuenta que se pueden simplificar un poco mas:

Es decir, los sistemas de ecuaciones lineales se pueden escribir como la multiplicación entre una matriz (que tiene los factores de las ecuaciones lineales) y el vector de incógnitas, igualada al vector de datos. Es decir, en general, los sistemas lineales se pueden escribir de la forma siguiente


Donde A es la matriz de factores, x el vector de incógnitas y b el vector con los datos. A los vectores los marqué con negritas, pero en el futuro no lo voy a hacer, excepto para enfatizar o para evitar confusión entre escalares y vectores.

Resolviendo sistemas de ecuaciones lineales

Cuando trabajamos con escalares (con números) sabemos como despejar una incógnita: se aplican las operaciones inversas hasta que la incógnita queda despejada. Si tenemos una ecuación lineal, podemos despejar de la siguiente forma

Como ven, multiplicamos ambos lados de la igualdad por la inversa multiplicativa de a. Pero para resolver el sistema de ecuaciones (1) ¿Cuál es la inversa multiplicativa de una matriz? La respuesta es un poco larga y no la vamos a discutir hoy (pueden investigar en la internet sobre "matriz inversa"). Lo bueno es que GNU Octave sabe como calcular la inversa de una matriz, así que van a poder resolver sistemas lineales antes de aprender la teoría. Con suerte esto los va a motivar para estudiar y entender la teoría.

Lo que si es necesario saber es que una matriz, a diferencia de un escalar, no siempre tiene inversa multiplicativa. En otras palabras ¡Hay problemas lineales que no tienen solución!

Pregunta 1:
¿Podés construir un sistema de dos ecuaciones lineales que no tenga solución?

Resolviendo el problema calculando la inversa

Si estamos seguros que el problema tiene solución, entonces la inversa de la matriz de factores tiene que existir y en Octave la calculamos con la función inv(). En nuestro ejemplo

A  = [2 0.5; -1 1];
b  = [10; 5];
Ai = inv (A);
x  = Ai * b

x =


Pregunta 2:
Si ahora tratás de calcular la inversa del problema imposible de la Pregunta 1 ¿Qué obtenés?

Resolviendo el problema sin calcular la inversa

Otra forma de resolver el problema en Octave es utilizar el operador \ (división izquierda). Este es un operador inteligente y elige la mejor forma de resolver el problema, desde un punto de vista numérico (ya veremos que significa esto). En la mayoría de los casos este operador no calcula la inversa explícitamente, puesto que dicha operación es un poco problemática. El código para resolver el sistema (1) con este operador es:

A = [2 0.5; -1 1];
b = [10; 5];
x = A \ b

x =


Pregunta 3:
Si ahora tratás de resolver el problema imposible de la Pregunta 1 ¿Qué obtenés?

¿Se chocan o no se chocan? Problemas de encuentro

Advertencia: el siguiente problema puede traerte recuerdos del colegio secundario :D.

Problema 2:
Un móvil A sale de un punto con una rapidez de 60 Km/h y luego de 30 minutos parte otro B del mismo punto a una rapidez de 80 Km/h ¿Cuándo alcanzará el segundo móvil al primero y a qué distancia del punto de partida se encuentran?.

Primero planteamos las ecuaciones de cada móvil:

Donde el subíndice indica a que móvil nos referimos. Estas ecuaciones nos dicen donde estarán los móviles A y B luego de un tiempo t. Vamos a reescribir estas ecuaciones para poner el problema en forma matricial (e.d. usando matrices). Primero notemos que nuestro objetivo es conocer en que tiempo t, las posiciones son iguales, e.d. xA = x= x
Reemplazando en las ecuaciones de arriba y poniendo las incógnitas del mismo lado obtenemos:

Ahora la matriz de factores es evidente y el problema se puede escribir como

El código de Octave que resuelve el Problema 2 es:

vA  = 60; t0A = 0; 
vB  = 80; t0B = 0.5;   % 30 minutos son media hora
A   = [ vA -1; vB -1 ];
b   = [ vA*t0A; vB*t0B ];
sol = A \ b

sol =

Los móviles se encuentran 2 horas luego de la partida de A (t=2) y se encuentran a 120 kilómetros del punto de partida (x = 120). Para verificar que la solución es correcta podemos reemplazar en las ecuaciones iniciales (o multiplicar A*sol y verificar que obtenemos b).

¡A resolver problemas!
Si tienen preguntas participen en el foro de discusión.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Producto escalar

En la clase sobre operadores vimos que el operador asterisco * parece no funcionar con vectores. Hoy veremos que este operador realiza el famoso producto escalar entre vectores o el producto entre matrices.

Longitud de un vector (módulo o norma)

Para indicar la posición de un punto en el plano (digamos una bolita en una cancha) necesitamos dar al menos dos magnitudes. Por ejemplo, tomando un punto como referencia (una esquina de la cancha), podemos indicar las distancias al punto desde ese origen en dirección horizontal y en dirección vertical. Estos dos valores los podemos acomodar en un vector [x y].

¿Cuál es la distancia en linea recta desde el origen hasta la bolita? La respuesta viene de la mano del teorema de Pitágoras

Este valor es la longitud del vector (t.c.c. módulo del vector o norma del vector). Dentro de la raíz cuadrada tenemos el cuadrado del módulo del vector, la suma del producto de cada componente por si misma. 
Con toda libertad podemos preguntarnos como se llama la misma operación entre dos vectores diferentes. Digamos que tenemos nuestro vector [x y] y un otro vector [u w]. Haciendo las cuentas obtenemos

Este valor es claramente un escalar (un número, no un vector), como lo es la distancia. A esta operación se la denomina producto escalar entre vectores. El producto escalar es una función que toma dos vectores como entrada y devuelve un valor escalar, de allí su nombre.

¿Matriz o vector?

Los matemáticos son tipos muy ocurrentes y se inventaron una forma de escribir el producto escalar de manera muy concisa. Primero decidieron diferenciar entre vectores fila y vectores columna (como ya vimos en otra clase) y luego decidieron que cuando uno multiplica un vector fila por un vector columna, se debe calcular el producto escalar. Es decir, multiplicamos el primer elemento del vector fila por el primer elemento del vector columna y lo sumamos al producto de los segundos elementos.

aimación producto escalar
Cambiar de vector filas a columnas en Octave es muy sencillo.

x  = [1 2];
xT = x'

xT =

y la operación se llama
transponer. Es decir la comillas simples ' es un operador unitario (toma una sola entrada) que devuelve el vector transpuesto: cambia filas por columnas o viceversa. Ahora podemos calcular el módulo cuadrado y la longitud del vector x definido más arriba:

mod2 = x * x';
long = sqrt (mod2)

Pregunta 1:
¿Qué hace la función sqrt()? ¿Por qué se llama así?

Ya nos vamos dando cuenta que no hay mucha diferencia entre un vector y una matriz. En efecto, el producto entre matrices es el producto escalar de los vectores fila de la primera matriz, por los vectores columna de la segunda.

A = [1 2; 1 -1];
# vectores filas
a1 = A(1,:); # primera fila
a2 = A(2,:); # segunda fila

B = [1 0; 0 1];
# vectores columnas
b1 = B(:,1); # primera columna
b2 = B(:,2); # segunda columna

C  = A * B;
Cp = [a1*b1 a1*b2; a2*b1 a2*b2];

assert (Cp == C)

La matriz Cp se construye usando el producto escalar de las filas y las columnas. La función assert produce un error si la condición que se pasa por argumento no es verdadera. Al no recibir ningún error, entendemos que las dos matrices son iguales.

¡Un vector es tan solo una matriz con una sola fila o columna!

¿Preguntas? Participa en el foro.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¿Verdadero o falso?

En el cole te enseñaron una cosa que llaman lógica formal (o solamente lógica). Si te gustó tenés suerte, probablemente tu profesor/a era muy buen@ y entendía el tema (¡yo tuve mucha suerte!). Lamentablemente la mayoría de los estudiantes con los que charlo odian la lógica. En esta clase voy a tratar de darles un ejemplo de porqué la lógica es super importante y que no es tan aburrida como muchos creen.

Esta clase puede ayudarte con los desafíos del nivel 1, si todavía no los resolviste.

"La lógica no sirve para nada"

¡Cuantas veces habrás escuchado esa frase o alguna otra similar! Como la mayoría de las herramientas matemáticas, la lógica no es más que una formalización de algo que hacemos naturalmente. Voy a tomar un ejemplo mundano: seleccionar fruta o verdura cuando vamos de compras.

Hace unos días, en uno de mis viajes de aprovisionamiento al mercadito de la esquina, observaba con cierta admiración a una mujer mayor seleccionar las manzanas que ponía en su bolsa de las compras. Su meticulosidad era sorprendente y me preguntaba si esa mujer no habría sido una matemática en sus años de actividad laboral.
Cada vez que la señora juntaba una manzana, la miraba de muchos ángulos diferentes, la golpeaba suavemente con un dedo, la olía y finalmente decidía si la manzana volvía a la góndola o si se convertía en parte de su compra.
Podemos hacer un modelito sencillo del proceso de selección y pensar que la señora ejecutaba una función que tomaba como entrada los aspectos de la manzana bajo observación. Esta función evalúa tres condiciones y decide si la manzana se compra o se devuelve:
  1. Tiene buena forma?
  2. Suena bien?
  3. Huele bien?
La respuesta a estas preguntas es si o no. Podemos re-escribir las preguntas como afirmaciones (o condiciones) y decidir si estas afirmaciones son verdaderas (en inglés: true) o falsa (en inglés: false):
  1. La manzana tiene buena forma.
  2. La manzana suena bien al ser golpeada suavemente.
  3. La manzana huele bien.
Si la manzana en cuestión hace que todas estas afirmaciones sean verdaderas, entonces la compramos. De lo contrario, si cualquiera de estas afirmaciones es falsa, devolvemos la manzana a la góndola.

Símbolos y operadores

En el ejemplo introduje los dos símbolos fundamentales de la lógica: true y false. Estos son los nombres que se utilizan en GNU Octave para estos valores, pero vos podés usar cualquier par de valores que te guste: (V,F); (0,1); (-1,1); (0,5); (blanco, negro); etc.. Para poder ejercer la lógica necesitamos saber como operar con estos valores.

Tomemos dos afirmaciones, una falsa y otra verdadera, para ejemplificar. 
  1. La afirmación verdadera es: la palabra "hacha" tiene 5 letras
  2. La afirmación falsa:  La palabra "hacha" no tiene ninguna letra "h".
Verifiquemos los valores lógicos de estas afirmaciones en GNU Octave. El siguiente código guarda el string "hacha" en una variable llamada p (de palabra) y evalúa la 1ra afirmación.

p = "hacha";
# Cuenta la cantidad de letras
nLetras = length (p);
# Evalua la 1ra afirmacion 
nLetras == 5

La primera linea de código guarda la palabra en la variable. La linea luego del primer comentario devuelve la longitud del string p, es decir la cantidad de letras, que guardamos en a variable nLetras (recuerda que las lineas que empiezan con el caracter # o %, son ignoradas por GNU Octave, son lineas de comentarios, y tu también puedes ignorarlas. Son solo ayuda). La última linea de código pregunta si nLetras es igual a 5. Notar que el operador == no asigna el valor a nLetras, sino que compara los valores y responde verdadero si estos son iguales (experimenta un poco con este operador!).
Si ejecutas esas líneas de código en una sesión de Octave obtendrás la respuesta 

ans = 1

En Octave el valor true se representa con un 1 (ejecuta true en tu sesión de Octave para verificar) por lo tanto el valor de la 1ra afirmación es verdadero. Veamos el código de la 2da afirmación:

all (p != "h")

¿Simple, eh? Dejo como ejercicio entender que esta ocurriendo aquí. Lo que necesitas saber es que all en inglés significa "todo/todos" y que el operador != es lo opuesto a == (verifica desigualdad en vez de igualdad) ¿Qué ocurre cuando comparamos el string p con el caracter "h"?

La respuesta de Octave a esta linea de código es 0. Indicando que, efectivamente, la segunda afirmación es falsa.

Ok, tenemos una afirmación verdadera y una falsa, y vamos a guardar estos valores de verdad en un vector:

afirmacion = [true, false]

De esta forma afirmacion(1) es verdadera (true, 1) y afirmacion(2) es falsa (false, 0).

¿Qué dirías de la nueva afirmación más abajo?
  • La palabra "hacha" tiene 5 letras y la palabra "hacha" no tiene ninguna letra "h"
¿Es verdadera o falsa?

El operador que hemos usado para construir esta nueva afirmación es el "y", en inglés "and", y que en lógica llamán conjunción. En Octave podemos utilizar la función and o el operador &, veamos: 

# operador
afirmacion(1) & afirmacion(2)
# funcion
and (afirmacion(1), afirmacion(2))

Ejecuta este código ¿Estas de acuerdo con Octave sobre el valor de verdad de la nueva afirmación?
¿Cuál sería el resultado si las dos afirmaciones fuesen verdaderas?

Veamos otra afirmación
  • La palabra "hacha" tiene 5 letras o la palabra "hacha" no tiene ninguna letra "h"
Para construir esta nueva afirmación utilizamos la disyunción inclusiva, el "o". Pero no el "o" exclusivo (que algunos escriben "ó" para que quede claro), sino un "o" que está contento si esto es verdad, si lo otro, o si ambos son verdad. Este operador en Octave se escribe | y la función es or:

# operador
afirmacion(1) | afirmacion(2)
# funcion
or (afirmacion(1), afirmacion(2))

Experimenta con este operador para entender como funciona. 

Los operadores & y |, al igual que las respectivas funciones toman como entrada dos afirmaciones, son operadores binarios (en el sentido que toman 2 entradas). Otros operadores binarios que ya conocés son el + y el *, ambos toman dos números como entrada y devuelven el resultado.

Ejercicio 1: ¿Puedes construir el resultado de los operadores & y | para todas las entradas posibles?

Ejercicio 2: ¿Puedes escribir la afirmación que la señora de la historia evaluaba para seleccionar las manzanas?

Ejemplo numérico

Ok, toda esta cháchara sobre lógica, pero ¿Cómo nos puede ser útil para los desafíos?

En los desafíos de nivel 1 es necesario filtrar o seleccionar ciertos números según sus propiedades. Por ejemplo ¿Cómo podemos seleccionar números pares?
Lo que sabemos es que un número par se escribe como

K = 2*N

donde N puede ser cualquier numero entero. Todo bien con esto, pero para cada número K que me dan, tendría que buscar un N, entre todos los enteros, que cumpla con esta condición...posible, pero medio pesado. De manera equivalente podemos decir que un número es par si al dividirlo en 2 el resto es cero. Es decir, un número par es divisible por 2. GNU Octave tiene una función que nos devuelve el resto de una división entera, la función rem (del inglés "remainder", que quiere decir "resto" o "sobrante") . Como ejemplo de uso tomemos la division de 7 en 2 partes. 7 se puede repartir entre 2 en 3 partes iguales, pero nos sobra 1. Veamos que dice a la función:

rem (7, 2)
ans = 1

El primer argumento es el numero a dividir (el dividendo) y el segundo argumento es el divisor.

Ahora vamos a utilizar esta función y nuestro conocimiento de lógica para seleccionar los números pares de una lista de números:

l = 0:3:15;
tf = rem (l,2) == 0
tf =
   1   0   1   0   1   0

La lista l contiene números enteros desde 0 hasta 15 en incrementos de a 3 y Octave nos dice que la afirmación "el número es divisble por 2" es verdadera para el primer número, falsa para el segundo, verdadera para el tercero y así sucesivamente. 
¿Cuales son estos números? Claramente 0, 6 y 12. En Octave podemos obtener estos números utilizando el vector lógico al que hemos llamado tf. Prueba lo siguiente


Sigue experimentando con este tipo de ejercicios hasta que entiendas como funciona la cuestión. Si tienes preguntas o sugerencias no dudes en escribirlas en el foro de discusión o los comentarios.

Adjunto un archivo con todos los comandos que hemos utilizado en esta clase.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Si, si en cambio, sino: IF-ELSEIF-ELSE

En la clase sobre lógica discutimos un poco la forma de generar valores verdadero/false, usualmente llamados valores booleanos (se lee bu-leanos y viene del nombre de un inglés llamado George Boole).

Casi todos los lenguajes de programación ofrecen una forma de controlar la ejecución de código utilizando como criterio valores booleanos. Es decir, cierta parte del código se ejecuta si una condición lógica es verdadera. Una forma ejemplar de realizar esta ejecución selectiva es utilizando IF-ELSE (si, sino) o IF-ELSEIF (si, si en cambio) o una combinación de ambos.


En Octave, todos aquellos comandos dentro de un bloque if se ejecutan solo si a condición dada es verdadera. 

if (x == 0)
  disp ("Equis es igual a cero");

El bloque if se empieza en la linea donde llamamos a esta función y termina con endif (o simplemente end, es lo mismo). En este caso, se imprime a la pantalla la oración "Equis es igual a cero" cuando la comparación x == 0 es verdadera. La condición puede ser tan compleja como sea necesario (aunque esto puede hacer el código difícil de entender), por ejemplo

if (x == 0 || x == 5)
  disp ("Equis es igual a cero o igual a cinco");

Si la condición fuese falsa, el código dentro del bloque nunca se ejecuta. En cualquier caso, una vez evaluado el bloque (con o sin ejecución del código correspondiente), el programa continua. El siguiente ejemplo imprime "Equis es igual a cinco" dependiendo del valor de x, pero siempre imprime "Terminé".

if (x == 5)
  disp ("Equis es igual a cinco");
disp ("Terminé.");

Pregunta 1:
¿En que situación se imprime a pantalla el mensaje "¡Muy loco!"?

if (x == 0 && x == 5)
  disp ("¡Muy loco!");

¿Veremos alguna vez el mensaje en pantalla?


El bloque if se puede extender para que cierto código se ejecute solo cuando la condición es falsa. Considera el siguiente código

disp ("Eqis es: ")

if (x > 0)
    disp ("positivo.");

if (x <= 0)
    disp ("negativo o cero.");

Dado que la segunda condición es lo opuesto a la primera podemos simplificar usando else

disp ("Eqis es: ")
if (x > 0)
    disp ("positivo.");
    disp ("negativo o cero.");

Pregunta 2:
¿Cómo harías para separar en "positivo", "cero" o "negativo" el ejemplo anterior?


Una respuesta a la pregunta anterior viene de la mano de elseif (si en cambio), que nos permite evaluar una condición adicional.

disp ("Eqis es: ")
if (x > 0)
    disp ("positivo.");
elseif (x == 0)
    disp ("cero.");
    disp ("negativo.");

Pregunta 3:
¿Podés dar otros ejemplos que produzcan lo mismo que el ejemplo anterior?

Debe notarse que el código dentro de un elseif se ejecutara solo y solo si la condición del if es falsa. En el siguiente ejemplo jamás veremos el mensaje "¡Ahá!" en pantalla aunque la condición del elseif es siempre verdadera.

if (true)
 disp ("Sí, es verdad.");
elseif (true)
 disp ("¡Ahá!");

Esto quiere decir que las condiciones en un bloque IF-ELSEIF deberían que ser mutuamente exclusivas (nunca ambas pueden ser verdaderas al mismo tiempo). Esto no es necesario, como se muestra en el siguiente ejemplo, pero la lógica resultante es confusa

if (x > 0 && x <=1)
    disp ("Entre cero y uno (incluido).");
elseif (x > 0.5)
    disp ("Más que un medio.");

Pregunta 4:
¿Podés simplificar la lógica del ejemplo anterior?

Espero recibir sus preguntas y dudas en el foro!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Tic,tac, tic, toc! ¿Cuánto tarda mi programa?

Matías Fajardo pidió que discutamos control de flujo de programas en GNU Octave. Es un pedido razonable así que vamos a ir en esa dirección. 
Con el tema de control de flujo tocamos algo importante para lenguajes como GNU Octave. Para entender la importancia del tema necesitamos primero poder medir el tiempo que tarde en correr un programa. Por esta razón vamos a empezar nuestro camino hacia control de flujo con unas funciones muy sencillas: tic y toc.


El objetivo es el siguiente: queremos medir cuanto tarda la ejecución de un programa. Por ejemplo ¿Cuánto tarda la creación de una matriz en Octave? Es decir ¿Cuánto tarda la ejecución del comando
rand (n)
para distintos valores de n?

Una forma de medir el tiempo de ejecución es usar las funciones tic y toc:
n = 100;
rand (n);
La función tic inicia un cronómetro y la función toc muestra le valor actual del cronómetro. 
En mi computadora los resultados para n igual a 1e2, 1e3, y 1e4 son
Elapsed time is 0.000303984 seconds.
Elapsed time is 0.024765 seconds.
Elapsed time is 1.57575 seconds.

Pregunta 1:
¿Qué tiempos obtenés vos?

La respuesta de toc nos dice el tiempo segundos que transcurrió desde la última vez que llamamos a la función tic. Este tiempo se llama tiempo de reloj de pared (wall-clock time, en Inglés) y corresponde más o menos al tiempo real que transcurrió en tu vida mientras se ejecutaba el código.

Pregunta 2:
¿Qué resultados obtenés si corres el siguiente código?
n = 5e2;
rand (n);

n = 5e2;
clear all
rand (n);
¿Notás algo raro?

¡Preguntas y comentarios en el foro!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

¡No soy cirujano pero necesito operar!

Así es, necesitamos operar, pero no en el sentido médico. Necesitamos operar con números.
En el ámbito de la programación llamamos operadores a un grupo de funciones que tienen una sintaxis particular (es decir que se escriben de una forma especial). Estas funciones generalmente coinciden con las bien conocidas operaciones de la aritmética:

 Operación Operador en Octave 
Suma  +
Resta  -
Multiplicación  *
División  / 

Generalmente, cuando queremos escribir "Calcule la suma entre el número 3 y 4" anotamos 3 + 4 (3 más 4). Pero podríamos escribir +(3,4) ("sume 3 y 4") ó + 3 4 (eso se llama notación polaca) o incluso 3 4 + ("3 y 4, ahora sume", notación polaca inversa). Lo importante es que la operación requiere dos entradas, el 3 y el 4, y el símbolo que representa la operación es el operador. Hay operadores binarios, como el que acabamos de ver, pero también los hay unitarios (solo una entrada) y ternarios (tres entradas)... y así sucesivamente.
La notación preferida, por ejemplo 3 + 4 para la suma, se llama notación de infijo (como prefijo, y sufijo) y el operador va entre medio de las entradas. Se darán cuenta que esta notación es muy clara para operadores binarios (con dos entradas).

En la teoría de hoy discutiremos otros operadores ademas de los aritméticos, que en Octave pueden resultar muy útiles: especilamente cuando queremos hacer operaciones entre vectores o listas de números. No voy explicar todos los operadores! Les explicaré unos cuantos, al resto lo pueden investigar y aprender por su cuenta.


Sumar dos vectores o listas de números

Supongamos que hemos colectado la duración, en minutos, del 1er tiempo de varios eventos deportivos en un vector fila llamado primerT

primerT = [20 45 30 30 46 22];

En otro vector hemos colectado la duración del 2do tiempo de los mismos eventos

segundoT = [27 49 35 30 45 20];

Y ahora nos gustaría saber la duración total del partido. Tenemos que sumar la duración de los dos tiempos. en Octave hacemos

TotalT = primerT + segundoT

Como pueden ver el operador suma se aplica elemento a elemento.

Multiplicar dos listas

Supongamos que en otra ocasión hemos colectado los precios del pasaje de colectivo de diferentes años

precio = [0.7, 1.2, 2.5, 4, 10];

Y en otro vector hemos colectado la cantidad de veces que viajamos en colectivo y pagado el precio correspondiente

veces = [10 8 5 5 2];

Ahora queremos calcular la plata que gastamos en cada situación. Tenemos que multiplicar los dos vectores. En Octave, si ingresamos veces * precio recibimos un error. 

Antes de seguir leyendo ¿Podés explicar el porqué?

El problema es que Octave considera a las dos entradas vectores, y entre vectores hay una operación llamada producto escalar (tcc. producto punto, producto interno) que en Octave se calcula usando el operador * (asterisco). Para indicarle a Octave que lo que queremos es un producto elemento a elemento utilizamos el operador .* :

veces .* precio

y obtenemos el resultado deseado.


Para elevar un número a una potencia, por ejemplo 3 al cuadrado, escribimos

3 ^ 2

el operador ^ representa la operación.

Para calcular el cuadrado de una lista de números operamos con .^


En general todo operador cuyo nombre empiece con un punto (.) realiza la operación elemento a elemento.

Pregunta 1:
¿Qué operador utilizarías para dividir dos listas de números?


Volvamos a las listas con duración de los tiempos de un deporte. Ahora vamos a organizar la información de otra manera que quizás es más fácil de visualizar.
Primero ordenemos los deportes asociando a cada uno con un número:
  1. Football
  2. Handball
  3. Basketball
Por cada deporte vamos a guardar la duración del 1er tiempo de diferentes eventos. A esta información la organizamos en forma de tabla:

   Football Handball  Basketball 
 1er Evento  45  20 15 
 2do Evento  49  20 18 
 3er Evento  46  25 18 
 4to Evento  45  22 15 

Podemos ingresar esta tabla en Octave de la siguiente manera:

primerT = [45 20 15; 49 20 18; 46 25 18; 45 22 15];

Supongamos que hacemos lo mismo para el segundo tiempo

segundoT = [49 21 15; 45 20 15; 46 20 15; 49 21 18];

La duración total la seguimos calculando con primerT + segundoT. Noten que esta operación se ejecuta, nuevamente, elemento a elemento.

Pregunta 2:
¿Cómo harías la multiplicación de estas dos tablas?

Hay muchísimos más operadores para investigar ¡Disfruten!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Archivando el trabajo

En esta clase veremos como crear y ejecutar un archivo m (mfile).

Todos los comandos que ejecutamos en GNU Octave pueden guardarse en un archivo de texto. Es importante que el archivo sea de texto plano (ascii), puesto que este es el formato que GNU Octave espera en los archivos de comandos. Esto es verdad en general para la mayoría de los lenguajes de programación.

Para empezar abrimos un archivo de texto con nuestro editor preferido (en mi caso, que trabajo en linux, uso Gedit). Enfatizo que no deben usar procesadores de texto como Office porque estos programas no guardan sus archivos en texto plano. Si estas en Windows puedes usar Notepad++.

El siguiente video muestra como trabajo en Ubuntu:

Crear y ejecutar mfile en octave-cli

El video muestra la creación de un mfile que corresponde a una función. Todavía no hemos aprendido a escribir funciones así que por el momento vamos a ver otro tipo de mfile: el archivo de comandos (script).

Abrí un archivo de texto y escribí los siguientes comandos:

disp ("Hola, este es un archivo de comandos! (script)");
a = 1:5;
txt = repmat (", %d",1,length(a)-1);
printf (["Los primeros %d números son: %d" txt "\n"],length(a), a);

Guarda este archivo notando el directorio donde quedará guardado (en el video utilicé el escritorio). Voy a suponer que el nombre que le diste al archivos es script_00.m
Ahora inicia una sesión de octave. Asegúrate que el archivo esta en el directorio actual: ejecuta pwd para ver el directorio actual y ls parta ver la lista de archivos.
Para ejecutar el script, simplemente lo llamamos por su nombre (sin extensión): 


¡Listo! Has creado tu primer script en Octave. La siguiente lista enumera los pasos para trabajar con scripts:
  • Abrir un archivo ASCII (texto plano) donde escribir los comandos.
  • Guardar el archivo con extensión .m (fijarse donde se guarda).
  • Llamar al archivo desde GNU Octave utilizando el nombre sin la extensión (asegurarse que pwd corresponde con el lugar done se guardó el archivo).
Es importante utilizar scripts (es la forma recomendada de trabajar), puesto que los comandos quedan archivados y podemos corregir errores y volver a ejecutarlos sin tener que escribir todo de nuevo.

Si estas usando una versión de Octave 3.8 o más nueva puedes utilizar la interfaz gráfica que incluye un editor de textos. El siguiente video muestra como crear un script que imprime los números del 1 al 10.

Crear y ejecutar mfile en octave-gui

La GUI de Octave puede configurarse en español si uds. lo desean. noten que esto solo cambia el idioma del menú, no el idioma del lenguaje de programación (funciones, ayuda y mensajes).

No duden en escribir sus preguntas y comentarios, aquí abajo o en el foro de discusión.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Crea tus propias funciones

A esta altura están en posición de escribir sus propias funciones. Pero ¿Qué es una función?

Una función es una entidad que recibe una serie de argumentos como entradas, opera sobre ellos y luego devuelve ciertos resultados. Las funciones cumplen un rol fundamental: encapsulan tareas. Al encapsular tareas, en especial aquellas que se repiten seguido, nuestros programas son mas fáciles de leer y a veces hasta más eficientes.

Considera la situación de saludar a personas que encuentras durante tu día. Supongamos que tu frase elegida como saludo es "Hola, <nombre>",  donde <nombre> es el nombre de pila de la persona a al que saludamos. Dado que a esto lo vamos a hacer muy seguido, podríamos definir una función (o rutina) que haga lo siguiente:

1. Toma el nombre de la persona a saludar.
2. Adjunta "Hola, ".
3. Imprima el saludo
4. Devuelva un resultado con el saludo.

Esta función en octave sería

function frase = saludar (nombre)
  frase = ["Hola, " nombre];
  disp (frase);

En una sesión de octave utilizaríamos esta función de la siguiente manera

x = "María";
y = saludar (x);

En la pantalla veríamos "Hola, María" y el contenido de la variable y sería exactamente esa frase.

Pregunta 1
¿Puedes definir una función que tome dos entradas y devuelva la suma de las mismas?

Estructura de una función

El siguiente esquema muestra la estructura de una función

function <resultados> = <nombre de función> ( <entradas> )

Las palabras y símbolos en negrita son necesarios y obligatorios.

Las palabra clave function es necesaria (y obligatoria) para definir una función. Un archivo que contiene una función debe siempre tener como primer comando ejecutable esta palabra clave. El archivo donde guardamos la función siempre tiene que llamarse  <nombre de funcion>.m , es decir el nombre del archivo debe coincidir con le nombre de a función que allí se define.
Las variables que se pasarán al ámbito desde donde se ha llamado la función (p. ej. la sesión de Octave) se devuelven como un vector en <resultados>. En el caso de la función saludar.m de más arriba, devolvíamos solo una variable. Pero una función puede devolver un montón de cosas o nada. Ejemplos

function [x, y] = ceroyuno ()
 x = 0;
 y = 1;

function diceHola ()
  disp ("Hola");

De manera similar, las variables de entrada  se pasan en una lista en <entradas>. Las entradas pueden ser muchas o ninguna (como en los ejemplos anteriores). Ejemplo

function z = resta (x,y)
  z = x - y;

Pregunta 2:
¿Cual es el nombre del archivo .m donde debemos guardar la función resta definida en el ejemplo anterior?

Pregunta 3:
¿Como se define una función que devuelva la cantidad de letras "a" en un string?


Cuando discutimos como guardar nuestro trabajo, vimos que desde Octave podemos ejecutar el contenido de un archivo de texto. Los comandos en ese archivo tienen acceso a las variables que existen en el ámbito (en inglés scope) de la sesión de Octave; esto quiere decir que los scripts pueden "ver" y "editar" las variables que creamos directamente en la sesión, y que si estos scripts crean nuevas variables estas quedaran en nuestra sesión luego de que el script haya terminado de ejecutarse.

Una función, en cambio, crea su propio ámbito de variables. Las variables que son creadas dentro de una función no son accesibles desde la sesión de Octave, o por lo menos no directamente. Excepto en casos particulares, las variables creadas dentro de una función son eliminadas cuando la función termina.

Con estas ayuditas ya puedes empezar crear tus propias funciones. El mundo de las funciones es mucho más amplio de lo que he descrito aquí ¡No te olvides de explorar!. Espero tus dudas y preguntas en el foro.

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

El que busca encuentra

El foro de discusión tiene una nueva entrada: una pregunta de Florencia que está resolviendo el desafío nro. 1.5.

La pregunta de Floppy nos da pié para discutir sobre un tipo de funciones muy útiles. Funciones que nos permiten obtener la posición, dentro de un vector, de elementos que cumplen cierta condición lógica. Esto es muy similar a lo que hicimos en la clase ¿Verdadero o falso? cuando buscamos números pares.

Find == Econtrar

Las operaciones lógicas generan matrices con valores verdadero o falso. Así, la siguiente desigualdad 

tf = [-1,5,-3] > 0

Genera un vector de valores lógicos (verdadero o falso), en este caso [false,true,false]. Estos vectores lógicos son de muchísima utilidad: nos permiten filtrar matrices o buscar elementos que cumplan ciertas condiciones lógicas. Para realizar lo segundo usamos la función find.
La función find (encontrar, en inglés) nos devuelve la posición de aquellos elementos que corresponden al valor "verdadero". Por ejemplo:

find ([false, true, false])

devuelve un 2, porque el segundo elemento es un valor verdadero. En Octave todo aquello que no es cero o vacío se considera verdadero. Es decir que

find ([0,-3, 0])

también devuelve 2, porque el segundo elemento es el único no cero. Cuando varios elementos son verdaderos, la función find devuelve todas las posiciones correspondientes, es decir

find ([false, true, true])

devuelve [2, 3].
La función también acepta matrices como entradas y nos devuelve el índice lineal de los elementos verdaderos. En Octave las matrices están ordenadas por columnas (en Inglés se dice column-major order), es decir contamos de arriba para abajo y de izquierda a derecha. Por ejemplo en una matriz de 3x4 los indices lineales son

Es decir que si ejecutamos

find ([0 1 0; 0 0 1; 0 1 0])

obtenemos [4,6,8]
Para obtener las filas y columnas de los elementos verdaderos, tenemos que llamar a la función con dos argumentos de salida. Si ejecutamos

[i, j] = find ([0 1 0; 0 0 1; 0 1 0]);
ans =
   1   2
   3   2
   2   3

obtenemos en i las filas y en j las columnas de los valores verdaderos.

¿Es miembro?

Otra función muy útil es ismember (es miembro, en inglés). Esta función toma como entrada dos argumentos y nos dice si los elementos del primero están en el segundo. Solo voy a proveer un corto ejemplo y los invito a que lean la ayuda de la función (en Octave: help ismember) y pongan sus preguntas en el foro de discusión.
Supongamos que queremos saber si la letra "a" está presente en un string

x = "Hay letras a en esta frase?";

podemos ejecutar

tf = ismember ("a", x)
  tf =  1

Claro que sí. La función también se puede usar para obtener la posición de las letras buscadas. Para esto intercambiamos los argumentos de entrada y buscamos el string x en la letra "a".

pos = find (ismember (x,"a"))
  pos =
      2    9   12   20   24

Pregunta 1. ¿Puedes explicar el porqué de este truquillo?

La función ismember también puede darnos la posición de los elementos encontrados, pero cuando el elemento buscado se repite (como en nuestro caso), solo nos devuelve la última ocurrencia

[tf, pos] = ismember ("a", x)
  tf  = 1
  pos = 24

Esta funcionalidad es útil cuando no hay repeticiones y se ilustra mejor con el ejemplo en la ayuda de ismember

a = [3, 10, 1];
s = [0:9];
[tf, pos] = ismember (a, s)
tf = [1, 0, 1]
s_idx = [4, 0, 2]

Las funciones find e ismember son de mucha utilidad y tienen varios otros modos de uso. Lee la ayuda, explora y pregunta en el foro!

by Juan Pablo Carbajal ( at July 18, 2016 07:56 PM

Cristiano Dorigo

Pcg and Gmres

Hi all,

during these two weeks I tried to fix some problems in pcg and gmres.
You can find the codes at
in the folder "pcg_gmres".


Question on my previous post

In my previous post I had some problem in pcg with the check if the matrix was or not positive definite. Indeed the check done in the Octave code was not precise, especially if involved complex numbers. (To recall: in the algorithm is computed alpha = (r' * r) / (p' * (A*p)) and the control was: if alpha <= 0, break. The problem is that also if alpha was complex the algorithm must stop, moreover for Octave every complex number is greater than 0, so pcg didn't notice if alpha has an imaginary part).

After a (brief) discussion with the Octave community I changed the control on alpha. First of all it checks separately the numerator and the denominator of alpha (because if there is a preconditioner matrix then the definition of alpha "slightly" changes and also the numerator can be negative and/or complex), then the algorithm stops if

(abs (imag (num)) >= eps (real (num))*tol/eps)  or  real (num) <= 0

(the same control is made for den).

In this way the algorithm stops if really the numerator (or denominator) of alpha is negative or if it has an imaginary part "too big" (this is relative to the tolerance, because we thought that it is not the case to stop the algorithm if the imaginary part is, for example, 1e-16). If this happen then the FLAG output is setted as 4.

Unfortunately, alpha can be real and positive at every iteration also if the matrix is not positive definite, so there are cases when this check doesn't works, but there are nothing that we can do.

Other changes

I made other changes to make the code more compatible with Matlab. The most relevant are:
  • Making some tests with Matlab I noticed that the output approximation X is not the last computed, but the one with the minimum residual and the output variable IT gives the iteration which X was computed. The Matlab help is not very explicit about this. I noticed this fact because controlling its outputs, sometimes lenght(RESVEC) is greater than IT, where RESVEC is the vector of the residuals. So at first I thought that there was something wrong (it was strange that there are 10, 20 or 30 residuals more than iterations), then I noticed this "minimum" relation. Then I modified the pcg in such a way also the Octave one has this output.
  • I added the check about the stagnation. Indeed there is the possibility that the algorithm stagnates (in general when the (preconditioned) matrix A is not positive definite but the algorithm doesn't notice it). So I add the control:
    if (norm (x - x_pr) <= eps*norm (x)) , where x_pr is the approximation of the previous iteration, then the algorithm stagnates and the FLAG output is setted as 3.
  • I added a check to control if the preconditioner matrix is singular, and if this happens then FLAG = 2. To do this, I use a try-catch approach:



    . . .   (pcg computations)

      flag = 2;

    In this way the warning "Octave:singular-matrix" is setted as an error, and with the try-catch I control if this error appears.
    I used this approach because in Matlab and Octave there are different results when we try to solve M \ b, with M singular:
    In Matlab there is always a warning and the output has some NaN or Inf;
    I Octave the warning appears only the first time that we try to solve it and the output is the least squares minimization.
    Unfortunately, since in Octave this warning appears only the first time, if we use pcg for a second time (without quitting Octave) with the same linear system with singular M, then FLAG = 2 is not setted.
    There is another case where this approach doesn't work: when the matrix is setted as "Diagonal Matrix" then this warning doesn't appear. For example:

    octave-cli:3> M1 = eye(3);
    octave-cli:4> M1(3,3) = 0
    M1 =

    Diagonal Matrix

       1   0   0
       0   1   0
       0   0   0

    octave-cli:5> M2 = [0 0 0;0 1 0; 0 0 1]
    M2 =

       0   0   0
       0   1   0
       0   0   1

    octave-cli:6> M1 \ ones(3,1)
    ans =


    octave-cli:7> M2 \ ones(3,1)
    warning: matrix singular to machine precision
    ans =


    In this case M1 and M2 are different singular matrices, but if I try to solve  a linear system with them there is a warning only with M2.
  • I setted the output flag as in Matlab:
    FLAG = 0 : algorithm converged
    FLAG = 1 : maximum number of iterations reached
    FLAG = 2 : preconditioner matrix is singular
    FLAG = 3 : stagnation
    FLAG = 4: (preconditioned) matrix A not positive definite

    About FLAG = 2 and FLAG = 4:

    "FLAG = 2": the Matlab help tell that this flag appears when M is ill-conditioned, but the only case in which I can reproduce this flag is when M is singular. But singular and ill-conditioned matrix are different concepts, (ill-conditioned when the condition number is high, singular when det(A) is zero). Indeed, for example, Hilbert matrices are ill-conditioned but not singular. Since in Matlab I had FLAG = 2 only when M singular, in the code and in my pcg help is written that this flag appears only when M is singular, and not when M ill-conditioned.

    "FLAG = 4": Also with this flag Matlab is not very clear. In the help is written that this flag appears when the algorithm can not continue due to some values too big or too small. But I can reproduce it only when A is not positive definite, so in my pcg help I wrote that it appears only when the algorithm detect that A is not positive definite.
  • I fixed the output strings that appears when there are asked less than 2 outputs and I wrote them in such a way they are similar to the Matlab ones.
  • I update the pcg help, in such a way it is synchronized with these changes.


More or less I made the same changes of the pcg algorithm also in gmres (the discussion about the positive definitiveness of the matrix doesn't care here, since the gmres algorithms must works with any matrix, so there are 3 types of flag and not 4 as in the pcg).
One little things about the variable IT. In both Matlab and Octave, there are written in the help that in the variable IT are memorized the inner and the outer iterations, but they are not so clear about what they are. For this reason in the help I added a little deep explanation about what they are and how they are related with the total number of iterations. In particular the inner iterations are the iteration before that the restart is applied, instead the outer iterations count how many times the restart is applied. Then the total number of iterations are:
(Outer Iterations)*RESTART + (Inner Iterations).

I think that the most of the work about these two algorithms is done. Up to small changes and/or bugs that I can find trying my codes I think that they are more or less definitive.

If you find some bugs, you find something not clear, or have some advices please contact me.

by Cristiano Dorigo ( at July 18, 2016 07:58 AM

July 14, 2016

Barbara Lócsi

chol/qz and matrix/vector flags

In the past weeks I worked on implementing the chol/qz and matrix/vector flags. My work can be found at my repository:

Calling forms

chol/qz flag

In case of a generalized eigenvalue problem the chol/qz flag can be used in Matlab. The aim of this task was to introduce this functionality to Octave.
The chol/qz flag only matters if your matrices are symmetric as:
Regardless of the algorithm you specify, the eig function always uses the QZ algorithm when A or B are not symmetric. [1]

If you omit the chol/qz flag it works as it worked before:
The algorithm used depends on whether there are one or two input matrices, if they are real or complex, and if they are symmetric (Hermitian if complex) or non-symmetric. [2]
Which is the same as in Matlab:
When you omit the algorithm argument, the eig function selects an algorithm based on the properties of A and B. It uses the 'chol' algorithm for symmetric (Hermitian) A and symmetric (Hermitian) positive definite B. Otherwise, it uses the 'qz' algorithm. [1]
So if it is a symmetric problem than the chol method is used otherwise the qz method is used.

matrix/vector flag

The matrix/vector flag can be used in every previous call forms. With it you can specify the format of the eigenvalues.
If you use the “vector” flag the eigenvalues will be returned in a column vector.
If you use the “matrix” flag the eigenvalues will be returned in a diagonal matrix.
If you omit this flag than the eigenvalues will be returned in the same format as before, which is the same as in Matlab:
The default behavior varies according to the number of outputs specified:
If you specify one output, such as e = eig(A), then the eigenvalues are returned as a column vector by default.
If you specify two or three outputs, such as [V,D] = eig(A), then the eigenvalues are returned as a diagonal matrix, D, by default. [3]


by Barbara Lócsi ( at July 14, 2016 02:54 PM

July 12, 2016

Chiara Segala

odexprk23 solver

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 described in my post. 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), in which there are three calls to expmv function. The first two for the calculation of the second and third stage
U_{n2} = u_n + c_2 h_n \varphi_1(c_2 h_n A) F(t_n,u_n) ,
U_{n3} = u_n + \frac{2}{3} h_n \varphi_1(\frac{2}{3} h_n A) F(t_n,u_n) + \frac{2}{3} h_n \varphi_2(\frac{2}{3} h_n A) \frac{2}{3 c_2}D_{n2} ,
myexpmv = @(t,A,v) expmv (t, A, v);

[Atilde, eta] = augmat (c2*dt, A, Fn);
X = myexpmv (c2*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:, 2) = x(:) + X;

[Atilde, eta] = augmat (2/3*dt, A, [Fn, 2/(3*c2)*D(:, 2)]);
X = myexpmv (2/3*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);

U(:,3) = x(:) + X;

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

u_{n+1} = u_n + h_n \varphi_1(h_n A) F(t_n,u_n) + h_n \varphi_2(h_n A) \frac{3}{2}D_{n3} ,
\hat{u}_{n+1} = u_n + h_n \varphi_1(h_n A) [ F(t_n,u_n) + \frac{1}{2 c_2}D_{n2}],
v = [Fn, 3/2*D(:, 3)];
w = [Fn + 1/(2*c2)*D(:, 2)];
[Atilde, eta] = augmat (dt, A, v, w);
p = size (v, 2);
q = size (w, 2);
ep = [[zeros(p-1, 1); 1/eta; zeros(q, 1)], [zeros(p+q-1, 1); 1/eta]];
X = [eye(d, d), zeros(d, p+q)]*myexpmv(dt, Atilde, [sparse(d, 1), sparse(d, 1); ep]);
X1 = X(:, 1);
X2 = X(:, 2);
x_next = x(:) + X1;
x_est = x(:) + X2;

I wrote also a code for the solver (odexprk23), according to odepkg standards.
I finally tested the order of the method with the usual example described in my post.

by Chiara Segala ( at July 12, 2016 02:03 AM

July 11, 2016

Amr Mohamed


This is my first post since the midterm evaluation so let’s celebrate it first.

During the past week , we improved the polybool script ,added comments and updated the documentation/help message .
Here is the most recent test that i have tried:

We also added two new .m files (poly2cw – poly2ccw) for adjusting the orientation of the polygons.
They are used to automatically convert polygons to clockwise/counterclockwise contours .

We are currently testing the polybool script and updating it to work smoothly on self-intersecting polygons.

We will also start working on the poly2fv function which will be the last required function that was planned for in the proposal.
Hopefully, We will manage to get things done so that we can work on more functions during the upcoming month.

by amrkeleg at July 11, 2016 08:38 AM

July 03, 2016

Cristiano Dorigo

Pcg analysis

Dear all,

In this second and half week I worked on the pcg to try to improve it.
All the codes can be found at   
 in the folder "pcg_codes".

The first thing that I made is the private function __checkandstring__ (the name is not definitive).
This function checks if the input argument A, b, M1 and M2 are consistent, i.e. if A, M1, M2 are matrices or function handle, and if A is a matrix, checks also that the number of columns of A are equal to the number of rows of b.
In addition, this function has as output two strings: is_prec and type.
The first one can be "prec" if M1 is not empty, or "no_prec" otherwise (I check only M1 because in the pcg documentation is written that if M1 is empty then no preconditioner is applied).
Instead the string type can be "matrix", "function_handle" or "mixed". If M1 empty then type is "matrix" if A is a matrix,  "function handle" otherwise.
If M1 is not empty then the string type is "matrix" if A, M1, M2 matrices, "function_handle" if they are all function handle or "mixed" if not all of them are matrices or function handle.
I made this function as general as possible to use it also for the other algorithms that need an improvement in my project.
You can see the code in the function __checkandstring__.m.

After that I take the pcg Octave code and I adapted it to the different situation:
I made a switch-case that treat differently the three cases mentioned before ("matrix", "function handle" and "mixed"). Moreover the cases "matrix" and "function handle", have an if block that check if there exist the preconditioners: if they are not passed, then applies the unpreconditioned conjugate gradient, otherwise applies the preconditioned one.
For the "mixed" case instead, differently from my previous post  (in which I said that if A, M1 and M2 are not of the same type, then I set them as functions handle) I check every time the type of A / M1 / M2 and if it is a matrix then I apply the simple matrix-vector product, if it is a function handle then I made the evaluation.
You can see this code in the function pcg_tmp_02.m.

I made this choice because I noticed that it is more (time) efficient than to transform all in functions handle. Indeed, before the two codes mentioned above, I made a function similar to __checkandstring__ that in case A, M1 and M2 not of the same type then transform them all in function handle, and then I made a version of pcg that works with this function (you can see these code in the files __check02__ and pcg_tmp_01 respectively). Then I tested the octave pcg, pcg_tmp_01 and pcg_tmp_02 20 times on the same example, with and without preconditioners, and passing all the variables as matrices, as function handles and not all of the same type. At the end I plotted the times and I compared it and I noticed that the better performs are given form pcg_tmp_02.
You can see this example in the file example_pcg.m.

Talking with the mentors, they told me that (using as example pcg_tmp_02) the code is too long. This because the pcg code is "repeated" 5 times:
  • matrix case without preconditioners
  • matrix case with preconditioners
  • function handle case without preconditioners
  • function handle case with preconditioners
  • mixed case (it is only with preconditiones because otherwise there is only A and it is or a matrix or a function handle)
They suggested me to substitute this switch-case approach with another approach, for example to use nested function or to insert some scripts in the private folder. They told me to ask to all the maintainers which is your preferred approach or to suggest me another one in case.

Complex case

Making some tests with complex matrices, I noticed strange behaviour of the pcg.
For example if I test the Octave pcg with the matrix

A = [2 , 2+1i,  4;
        2-1i , 3, 1i;
        4 , -1i , 1 ]

that is an hermitian matrix, but not positive definite (indeed its eigenvalue are not all positive:   -3.16512506440744e+00   2.85307928516935e+00  6.31204577923809e+00 )
and with right-hand-side

b = [ 5.59693687377178e-01;
         9.77795880620014e-01 ]

the Octave pcg doesn't notice that this matrix is not positive definite.
Indeed the result is:

x = pcg(A,b)
warning: pcg: maximum number of iterations (3) reached
warning: the initial residual norm was reduced 2.12736e+15 times.
x =

   2.04219909327648e-01 - 3.21781116721895e-02i
  -1.11630690554651e-01 + 5.38299446813016e-02i
   1.07086298628122e-01 + 1.70817561341068e-02i
Studying the code, I notice that the only check for the positive definitiveness of the matrix is on alpha. Not going into deep details, this alpha is necessary to compute the approximant solution at every step (indeed at every iteration, x = x + alpha*p for a certain p).
If there are no preconditioners (as in our example) alpha is defined as:
alpha = < r, r > / < p, A*p >
(where < , > is the  vector scalar product, and r is the residual of the previous iteration).
Then, since we suppose that A is positive definite, alpha must be real and positive (both in case A real or complex matrix).
Since the numerator is positive for every r, if alpha is negative then surely A is not positive definite. Indeed the check in the Octave pcg code is:

 if (alpha <= 0.0)
      ## Negative matrix.
      matrix_positive_definite = false;

If A is not positive definite then < p, A*p> is not always negative, so if we are "lucky", if we pass a not positive definite matrix in pcg and from the computation alpha is positive at every iteration, there is nothing that we can do.
But if we use A and b definite as previous, pcg makes three iterations and the alpha are:

alpha =  2.11061741184511e-01 - 6.08511684384459e-20i
alpha = -5.59991838158184e-01 + 1.05798887633281e-17i
alpha =  1.48434184396776e-01 + 8.13133738226767e-18i

(the function pcg01.m in the folder "pcg_codes" is the octave version of pcg that print at every iteration alpha)

We notice that these three alpha have an imaginary part, but it is under the machine precision, so they are not computationally relevant.
But the second alpha is (obviously) negative, and the pcg doesn't notice it!

This because for Octave:

octave4.1-cli:29> -5.59991838158184e-01 + 1.05798887633281e-17i < 0
ans = 0


octave4.1-cli:30> -5.59991838158184e-01 + 1.05798887633281e-17i > 0
ans =  1

Making some tests I notice that a number with an imaginary part is greater than 0, also if this imaginary part is negative, indeed

octave4.1-cli:31> -3-1i > 0
ans =  1

I don't understand very well what Octave do to compare two complex numbers (I think that it compares the modules, i.e. if we want to verify if a < b, where at least one of a or b is complex, then it compares abs(a) and abs(b), but I'm not sure).
I compared some complex numbers also in Matlab, and I think that it compares only the real part (and "forgetting" the imaginary part).

I talked with the mentors about this fact, but there are some questions:
  • Is it correct to check only the real part of alpha?
  • Also if the imaginary part of alpha is relevant?
    (if for example we use the matrix

    A = [ 0.71290 + 0.59353i   0.97470 + 0.36591i   0.50060 + 0.53652i;
       0.37411 + 0.11662i   0.38904 + 0.43489i   0.03555 + 0.23431i;
       0.35482 + 0.23601i   0.44859 + 0.31402i   0.54356 + 0.72676i]

    that is a matrix not symmetric, and not positive definite,
    and b = ones(3,1), we obtain as alpha:

    alpha =  0.47882 - 0.51052i
    alpha =  1.4276 - 1.0314i
    alpha =  0.43535 - 0.36971i

    and the pcg doesn't give any flag of this "non-correctness" of A, but also if we check only the real part of alpha the code doesn't give any flag or any error)
  • Do we break the algorithm if alpha is complex? Also if the imaginary part is under the machine precision?
    (it is possible that also if A is positive definite and symmetric but the alpha are complex, because of computations and the machine precision)
The mentors suggested me two things:
  • First of all to check both numerator and denominator of alpha: this because if there are preconditioners also the numerator can be negative, since it becomes <M*r, r>, where M is the inverse of M1*M2.
  • To break the algorithm if the imaginary part of (the numerator or the denominator of) alpha is relevant. But, how to decide if it is relevant?
    For example, 1e-15 is relevant? And 1e-14?
    So they suggested me a possible criterion to use to determine the relevancy of this imaginary part:

    if real(alpha) + (eps / tol)*imag(alpha) == real(alpha)

    The motivation of this criterion is because, since we want to find a solution within a tolerance (tol), we "normalize" the relevancy of this imaginary part with this tolerance.
    The mentors also told me to ask also to all the maintainers, because this is their suggestion, but can be different from the direction that wants the Octave community.
    Then I please ask from the community some advices of how to proceed for this situation.

    Actually (and until I'll have some advices) I'm working on the gmres, and I'll try to improve it in the same manner of the first part of this post for the pcg.

    Please contact me with some feedback via e-mail at or as a comment under this post.

    by Cristiano Dorigo ( at July 03, 2016 03:59 AM

    June 25, 2016

    Francesco Faccio

    Midterm week and next steps

    Hello! This is my final post before the end of midterm evaluation and I would like to show what I've done during these days and which are the next steps. After the previous post I wrote some of the options that the user can provide using odeset. I wrote a function to allow the user to pass the modified Jacobian through a function that returns DF/DY and DF/DYP or through a cell array (if the Jacobian is constant) and I completed options MaxStep, InitialStep, MaxOrder (commit ea1a311). I started with mentors to write function decic, which computes initial conditions of the problem F(t,x,x') considering some values of x and x' fixed by the user. This function will be submitted as a patch in the next weeks. I tried to build Octave with KLU module included (it is an optional part of IDA which depends on SuiteSparse) but I had some problems with CMake that I hope to solve during the next days. KLU module will be used to solve sparse problems in ode15i, when a sparse Jacobian is supplied. Once I will be able to write the "sparse part" of ode15i, I will show some efficiency tests between my .oct file, Sundials' MEX and C file and ode15i in Matlab. After completing this part, I will try to improve the quality of the code: - writing a class whose methods will be able to call the function supplied by the user without global pointers and perform more efficient data conversion (avoiding loops) - Refactoring the code, putting all the checks in an other function For the examples I provided in the previous post, consider commit 4f60e96.

    by Francesco Faccio ( at June 25, 2016 07:13 PM

    June 21, 2016

    Barbara Lócsi

    Midterm evaluations

    You can find my work here:

    Project goals

    Certain calling forms of the eig function are currently missing, including:
    • preliminary balancing
    • computing left eigenvectors as a third output
    • choosing among generalized eigenvalue algorithms
    • choosing among return value formats of the eigenvalues (vector or matrix) see more here.

    Calling forms

    The aim for this period was:

    • Finish preliminary balancing if it is not finished, start working on implementing left eigenvector calculation

    The preliminary balancing task is done, I wrote a blog post about it:
    I started working on implementing the left eigenvector calculation and it is mostly done, some tests needs to be added to which will be completed until the end of this week.

    Remaining tasks for the second period:

    1. choosing among generalized eigenvalue algorithms
    2. choosing among return value formats of the eigenvalues (vector or matrix) see more here.

    by Barbara Lócsi ( at June 21, 2016 02:33 AM

    June 20, 2016

    Chiara Segala

    Week 4: augmented matrix

    In this week I wrote a code that implements the augmented matrix 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. The code is

    function [Atilde, eta] = augmat(h,A,V)

    p = size(V,2);
    W = fliplr(V/diag(h.^(0:p-1)));
    eta = 2^-ceil(log2(max(norm(W,1),realmin)));
    Atilde = [ A , eta*W ; zeros(p,size(A)) , diag(ones(p-1,1),1) ];

    in [HAM 11] eta is introduced to avoid rounding errors. I made a small change to avoid further errors, I added the max function so that smaller values ​​than realmin will not be considered.

    Now I summarize my work during this month.

    Week 1: phi functions. I implemented four functions, based on [BSW 07] (phi1m.m, phi2m.m, phi3m.m, phi4m.m).

    Week 2-3: general exponential schemes. I implemented the two schemes for a general exponential Runge-Kutta and Rosenbrock integrator (exprk.m, exprb.m).
    As already mentioned, these schemes are not really fast and efficient, but these are working codes that I applied to four different exponential methods (exprk2.m, exprk3.m, exprb3.m, exprb4.m) and I will use them as a reference when I go to implement the official methods.

    Week 4: augmented matrix (augmat.m).

    My codes can be found here

    by Chiara Segala ( at June 20, 2016 12:02 PM

    Cristiano Dorigo

    First week - Studying codes and some problems

    Dear all,

    In this first and half week, as written in the first post, I studied the algorithms and also their codes in Octave. I recall that the algorithms are gmres, pcg, bicg, bicgstab, qmr, cgs and pcr.

    The codes that I used that are not in Octave can be found at
    in the folder "First_week".

    Octave codes

    Studying the implementations I noticed problems in some of them:

    1. pcg - Trying to solve the simple problem 1*x = 1, with initial guess 1, there is the following output:

      octave4.1-cli:29> x = pcg(1,1,[ ],[ ],[ ],[ ],1)
      warning: division by zero
      warning: called from
          pcg at line 367 column 10
      pcg: converged in 0 iterations. pcg: the initial residual norm was reduced NaN times.
      x =  1

      We can see that in 0 iteration there is a warning for a division by zero and that the residual is reduced NaN times. Then this output text needs to be fixed.
    2. bicg - I tried to solve the small linear system A*x = b, with
      A = [7 7 8 4 4; 7 0 6 3 4; 3 2 3 7 6; 6 0 9 7 7; 1 0 0 1 7];
      b = sum(A,2);
      (I obtained this matrix with floor(10*rand(5)) so it is a simple integer matrix, and b computed as the sum of the columns of A in such a way the exact solution is the vector x = [1; 1; 1; 1; 1])
      With this configuration I discovered that the algorithm stops before the convergence, because of stagnation. Indeed the output is:

      octave4.1-cli:39> [x, flag, relres, iter, resvec] = bicg(A,b)
      bicg stopped at iteration 4 without converging to the desired tolerance 1.000000e-06
      because the method stagnated.
      The iterate returned (number 3) has relative residual 3.772464e-02
      x =


      flag =  3
      relres =  0.038052
      iter =  4
      resvec =

         0.070581   0.042765   0.037725

      As we can see the last residual (relres = 0.038052) is greater than the third (that is the third component of resvec, i.e. 0.037725). Then the algorithm detect a stagnation and stops.
      I tried to solve the same linear system with Matlab, and I noticed that also its bicg has the fourth iteration with the residual greater than the third one, but it doesn't stop and at the fifth iteration it converges at the solution.
      So I tried to comment the part of the stagnation (i.e. lines from 187 to 195 in the Octave's bicg) and tried to solve this linear system, obtaining:

      octave4.1-cli:41> [x,flag,relres,iter,resvec] = bicg_1(A,b)
      bicg stopped at iteration 5 without converging to the desired tolerance 1.000000e-06
      because the maximum number of iterations was reached. The iterate returned (number 5) has relative residual 3.371209e-14
      x =


      flag =  1
      relres =    3.3712e-14
      iter =  5
      resvec =

         0.070581   0.042765   0.037725   0.038052

      Except that an error text that says that the maximum number of iterations was reached, we can see that at the fifth iteration the algorithm converges with tolerance smaller than the default one 1e-6.
      The bicg with the "stagnation" commented is the file bigc_1.m, instead the matrix A and the test with bigc and bicg_1, can be found in example_bicg.m, in the folder First_week .
    3. cgs - I studied the algorithm from the Saad's book, and then I tried to study the codes, but I noticed that its Octave implementation was not taken from this book. Since in the help there is not the reference book I can't compare this implementation with the reference.
      I also tried to solve the same linear system discussed in the bicg section, and I noticed that the algorithm does not converge. The obtained output is:

      octave4.1-cli:56> [x,flag,relres,iter,resvec] = cgs(A,b)
      x =


      flag =  1
      relres =  0.17700
      iter =  5
      resvec =


      We can see that the solution given by the algorithm is far from the exact solution, also the flag equal to 1 indicates that the algorithm reached the maximum number of iterations without convergence.
      Then I tried to make a simple and basic implementation from scratch using the formulation in the Saad book, I tried to solve the linear system with this implementation and this implementation converges. The output is:

      octave4.1-cli:58> [x,resvec] = cgs_1(A,b)
      x =


      resvec =


      To be sure that is not a pathological case I tried to solve with the Octave implementation and mine the linear system A*x = b, with
      A = [6 6 8 7 2; 5 4 1 0 3; 4 1 2 6 3; 6 9 3 8 3; 3 5 6 9 4];
      b = sum(A,2);
      obtaining similar results.

      Then, after a discussion with the mentors, we decided that maybe it is easier to implement this method from scratch, instead to search the problems, due to the fact that there is no reference for its implementation.
      My simple and basic cgs implementation is the file cgs_1.m, and the test with the second matrix can be found in the file example_cgs.m, in the folder First_week.
    4. bicgstab - I noticed that in Matlab are displayed also the "half iteration" with the "half residual", but not in Octave. Then talking with the mentors we decided to display them also in Octave.

     Matrix-vector product in the algorithms

    Since all of the algorithms of the project are solvers for linear systems, the matrix-vector product is essential in all of their implementation, but since the matrix A and the (eventual) preconditioners can be passed as matrices or as function handle, it is not easy to decide how to compute all these matrix-vector products.

    For example in the gmres algorithm the strategy is to transform A, M1, M2 into function handle, i.e. (tanking A as example) to create the function handle

    Ax = @(x) A*x;  if A is a matrix
    Ax = A;              if A is a function handle

    and then every time that there is a matrix-vector product, for example A*v, it is simply computed as Ax(v) or feval(Ax, v).

    Instead, in the pcg implementation, every time that there is a matrix-vector product it is checked if A is a matrix or a function handle. For example if we need to compute A*v then there are the following lines of codes:

    if (isnumeric(A))
      w = A*v;
      w = feval(A,v);

    (the other algorithms use one of these two strategies, or similar variations).

    Then it is natural the question: which is the best strategy?
    To answer, I simulate these strategy and the simple matrix-vector product, then  I computed the times, using sparse matrices of different dimensions (n = 100,
    n = 200, n = 500, n= 1000, n = 2000) in this way 

      A = sprand(n, n, rand(1));
      Ax = @(x) A*x;
      w = rand(length(A),1);

     # Estimation of the simply A*w
      for j =1:100
        u = A*w;
      t1 = toc;

     # Estimation of feval(Ax, w)
      for j = 1:100
        u = feval (Ax, w);
      t2 = toc;

    # Estimation of the if block with A as matrix
      for j = 1:100
       if (isnumeric (A))
         u = A*w;
         u = feval(A,w);
      t3 = toc;

    # Estimation of the if block with Ax the function handle
        for j = 1:100
       if (isnumeric (Ax))
         u = Ax*w;
         u = feval(Ax,w);
      t4 = toc;

    obtaining (for example, since I used random matrices) as results:

    n = 100
    t1 =  0.0038838
    t2 =  0.0068600
    t3 =  0.0049579
    t4 =  0.0076430

    n = 200
    t1 =  0.0051990
    t2 =  0.0083339
    t3 =  0.0062981
    t4 =  0.0091600

    n = 500
    t1 =  0.0042169
    t2 =  0.0071800
    t3 =  0.0051830
    t4 =  0.0081360

    n = 1000
    t1 =  0.055708
    t2 =  0.059194
    t3 =  0.056590
    t4 =  0.060054

    n = 2000
    t1 =  0.15554
    t2 =  0.16691
    t3 =  0.16327
    t4 =  0.16677

    This test with different matrices can be found in the file quicktest.m in the folder First_week.

    We can see that if we consider n<=1000, the best strategy is to use the simple matix-vector product (indeed t1 is the less value). Instead if we consider the two strategies of gmres and pcg, we notice that the "if block" used in pcg takes more time than the simple matrix-vector product and the evaluation of the function handle.
    Instead the case n = 2000 is not a good example to see the differences of these approaches, probably because the matrix is too big and requires a lot of operations.
    Then, after a discussion with the mentors, we decided that the "best strategy" is to write in every algorithm two nested function "ad hoc" to the type of the data, and if A, M1 and M2 are not of the same type, to transform all of them in function handle. To be clear, for example, as sketch of how it will be:

    function [x, flag, relres, iter, resvec] = pcg (A, b, tol, maxit, M1, M2, x0)

      if A, M1, M2 are matrices
        [x, flag, relres, iter, resvec] = pcg_matrix (A, b, tol, maxit, M1, M2, x0)

      if A, M1, M2 are function handle
       [x, flag, relres, iter, resvec] = pcg_function (A, b, tol, maxit, M1, M2, x0)

      if  A, M1, M2 are not of the same type
       [x, flag, relres, iter, resvec] = pcg_mixed (A, b, tol, maxit, M1, M2, x0)


     Priority of the algorithms

    The last part of this post is about the chosen priority of the improvements of the algorithms for the next weeks.
    The mentors suggested me to improve the algorithms in the following order:
    1. pcg
    2. gmres
    3. bicgstab
    4. pcr
    5. qmr
    6. cgs
    7. bicg
    This order is suggested because pcg, gmres and bicgstab are the most used to solve iteratively linear systems.
    Then follows pcr because it is the only method already implemented in Octave that solve linear system when the input matrix A is Hermitian (bicgstab and gmres work for general matrices, instead pcg works for SPD matrices).
    And in the last follows qmr, cgs and bicg because they are the less used methods, with cgs above bicg because it has to be written from scratch.

    If there are something that is not clear or if there are some suggestions/advices, comment below or contact me via email at

      by Cristiano Dorigo ( at June 20, 2016 07:45 AM

      Francesco Faccio

      Summary of the work of the first month


      Mid-term evaluation has now arrived so it's time to summarize the work I've done and check which goals I have achieved. During this period I enjoied working with the community and the advices given by the mentors and by the other members have been really helpful. The most important change to the project is that, discussing with mentors, we decided to start implementing ode15i because it's more general than ode15s and to build ode15s later around it.
      Here you can find the code I've written so far:

      The most difficult task in the first part of this project was to have Octave compiled with link to Sundials. After having accomplished this, I checked the presence and usability of the library nvector_serial which contains the implementation of IDADENSE and IDABAND modules. I aggregated its build flags with the flags of sundials_ida and included the header nvector_serial.h in the dld-function.

      I checked the licenses of Sundials and SUPERLUMT (a package which will be used as a sparse direct solver, independent from Sundials) and they have 3-Clause license, so they are compatible with GNU-license and can be used.

      After configuring the correct flags, I began writing a minimal wrapper of ode15i of the form:

      [t, y] = ode15i (odefun, tspan, y0, yp0, options)

      The first problem was to deal with Sundials' own types. Sundials uses realtype and N_Vector. An N_Vector is a vector of realtype, while a realtype can be both a float, a double or a long double, depending on how Sundials has been built (default type is double). I assumed to use the default double realtype and wrote functions N_Vector ColToNVec (ColumnVector data, long int n) and ColumnVector NVecToCol (N_Vector v, long int n), which convert an Octave ColumnVector to an N_Vector and viceversa.

      I checked some minimal input conditions, wrote a few input validation tests, set AbsTol, Reltol, tspan, y0 and yp0 of type realtype or N_Vector.

      Once preprocessed data, the moment of glory of Sundials' functions arrived. The first call was to IDACreate(), which creates an IDA memory block and returns a pointer which is then passed as the first argument to all subsequent ida function calls.

      Sundials then asks to provide a function which computes the residual function in the DAE. This function must have the form:

      int flag = resfun (realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user data)

      As a temporary workaround I wrote a function which converts yy, yp and tt in ColumnVector(s), uses feval to evaluate the DAE (passed through a global pointer of type octave_function) and put the output in rr.

      Then a call to IDAInit, IDASVtolerances (or IDASStolerances if AbsTol is a scalar), IDADense and IDADlsSetDenseJacFn (if supplied) sets up the linear solver.

      Sundials accepts only a Jacobian Function of the form  J = DF/DY + cj*DF/DYP, where cj is a scalar proportional to the inverse of the step size (cj is computed by IDA's solver). I used the same workaround of the residual to evaluate J when it's required.

      Finally in the main loop a call to IDASolve solves the DAE system and gives the solution in output.

      What this wrapper can solve:

      this wrapper of ode15i can solve a system of differential equations of the form f(t, y, y') = 0, integrating from t0 to tf with initial conditions y0 and yp0. The output of the function is the solution of the DAE evaluated ONLY in the points supplied in tspan.
      It accepts as option a scalar RelTol and a scalar or vector AbsTol. If the user wants to supply a Jacobian, it must be of the form J = DF/DY + cj DF/DYP. Both the system of DAE and the Jacobian must be a function handle.

      I tested this wrapper using the 2 benchmark problems described in the previous post.

      In Robertson chemical kinetics problem I found the right solution both passing a Jacobian and letting Sundials approximate it. The script I used is the following:

      function res = robertsidae(t, y, yp)
      res = [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3));
      -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2);
      y(1) + y(2) + y(3) - 1];

      function jacc = jacobian(t, y, yp, c)
      jacc = [-0.04-c, 1e4*y(3), 1e4*y(2);
      0.04, -c-6*1e7*y(2)-1e4*y(3), -1e4*y(2);
      1, 1, 1];

      options = odeset('RelTol', 1e-3, 'AbsTol', 1e-6, ...
      'Jacobian', @jacobian);

      y0 = [1; 0; 0];
      yp0 = [-1e-4; 1e-4; 0];
      tspan = [0 4*logspace(-6, 6)];

      [t, y] = ode15i(@robertsidae, tspan, y0, yp0, options);

      y(:, 2) = 1e4*y(:, 2);
      semilogx(t, y)
      ylabel('species concentration');
      title('Robertson DAE problem with a Conservation Law');
      legend ('y1', 'y2', 'y3');

      As a result I got this:

      The second problem was to find the solution of a 2-D heat equation semidiscretized to a DAE on the unit square, as described in the previous post.

      After discretizing the domain, the system of 100 differential algebraic equations was passed to ode15i (as sparse methods are in progress, I used the dense solver of Sundials also for this problem, without passing the Jacobian).

      Here you can find the script I wrote to solve the problem. A 3-D time dependent plot shows how the solution evolves in time and space.

      uu0 = zeros(100, 1);

      %Initialize uu in all grid points
      for j = 1:10
      yfact = (1 / 9).*(j - 1);
      offset = 10.*(j - 1);
      for i = 1:10
      xfact= (1 / 9).*(i - 1);
      loc = offset + (i - 1);
      uu0(loc + 1) = 16 * xfact * (1 - xfact) * yfact * (1 - yfact);

      up0 = zeros(100, 1);

      %Set values of uu and up at boundary points
      for j = 1:10
      offset = 10 * (j - 1);
      for i = 1:10
      loc = offset + (i - 1);
      if (j == 1 || j == 10 || i == 1 || i == 10 )
      uu0 (loc + 1) = 0; up0 (loc + 1) = 0;

      function res = klu (t, uu, up)
      res = [ uu(1);
      up(12) - 81 * (uu(11) + uu(13) + uu(2) + uu(22) - 4*uu(12));
      up(13) - 81 * (uu(12) + uu(14) + uu(3) + uu(23) - 4*uu(13));
      up(14) - 81 * (uu(13) + uu(15) + uu(4) + uu(24) - 4*uu(14));
      up(15) - 81 * (uu(14) + uu(16) + uu(5) + uu(25) - 4*uu(15));
      up(16) - 81 * (uu(15) + uu(17) + uu(6) + uu(26) - 4*uu(16));
      up(17) - 81 * (uu(16) + uu(18) + uu(7) + uu(27) - 4*uu(17));
      up(18) - 81 * (uu(17) + uu(19) + uu(8) + uu(28) - 4*uu(18));
      up(19) - 81 * (uu(18) + uu(20) + uu(9) + uu(29) - 4*uu(19));
      up(22) - 81 * (uu(21) + uu(23) + uu(12) + uu(32) - 4*uu(22));
      up(23) - 81 * (uu(22) + uu(24) + uu(13) + uu(33) - 4*uu(23));
      up(24) - 81 * (uu(23) + uu(25) + uu(14) + uu(34) - 4*uu(24));
      up(25) - 81 * (uu(24) + uu(26) + uu(15) + uu(35) - 4*uu(25));
      up(26) - 81 * (uu(25) + uu(27) + uu(16) + uu(36) - 4*uu(26));
      up(27) - 81 * (uu(26) + uu(28) + uu(17) + uu(37) - 4*uu(27));
      up(28) - 81 * (uu(27) + uu(29) + uu(18) + uu(38) - 4*uu(28));
      up(29) - 81 * (uu(28) + uu(30) + uu(19) + uu(39) - 4*uu(29));
      up(32) - 81 * (uu(31) + uu(33) + uu(22) + uu(42) - 4*uu(32));
      up(33) - 81 * (uu(32) + uu(34) + uu(23) + uu(43) - 4*uu(33));
      up(34) - 81 * (uu(33) + uu(35) + uu(24) + uu(44) - 4*uu(34));
      up(35) - 81 * (uu(34) + uu(36) + uu(25) + uu(45) - 4*uu(35));
      up(36) - 81 * (uu(35) + uu(37) + uu(26) + uu(46) - 4*uu(36));
      up(37) - 81 * (uu(36) + uu(38) + uu(27) + uu(47) - 4*uu(37));
      up(38) - 81 * (uu(37) + uu(39) + uu(28) + uu(48) - 4*uu(38));
      up(39) - 81 * (uu(38) + uu(40) + uu(29) + uu(49) - 4*uu(39));
      up(42) - 81 * (uu(41) + uu(43) + uu(32) + uu(52) - 4*uu(42));
      up(43) - 81 * (uu(42) + uu(44) + uu(33) + uu(53) - 4*uu(43));
      up(44) - 81 * (uu(43) + uu(45) + uu(34) + uu(54) - 4*uu(44));
      up(45) - 81 * (uu(44) + uu(46) + uu(35) + uu(55) - 4*uu(45));
      up(46) - 81 * (uu(45) + uu(47) + uu(36) + uu(56) - 4*uu(46));
      up(47) - 81 * (uu(46) + uu(48) + uu(37) + uu(57) - 4*uu(47));
      up(48) - 81 * (uu(47) + uu(49) + uu(38) + uu(58) - 4*uu(48));
      up(49) - 81 * (uu(48) + uu(50) + uu(39) + uu(59) - 4*uu(49));
      up(52) - 81 * (uu(51) + uu(53) + uu(42) + uu(62) - 4*uu(52));
      up(53) - 81 * (uu(52) + uu(54) + uu(43) + uu(63) - 4*uu(53));
      up(54) - 81 * (uu(53) + uu(55) + uu(44) + uu(64) - 4*uu(54));
      up(55) - 81 * (uu(54) + uu(56) + uu(45) + uu(65) - 4*uu(55));
      up(56) - 81 * (uu(55) + uu(57) + uu(46) + uu(66) - 4*uu(56));
      up(57) - 81 * (uu(56) + uu(58) + uu(47) + uu(67) - 4*uu(57));
      up(58) - 81 * (uu(57) + uu(59) + uu(48) + uu(68) - 4*uu(58));
      up(59) - 81 * (uu(58) + uu(50) + uu(49) + uu(69) - 4*uu(59));
      up(62) - 81 * (uu(61) + uu(63) + uu(52) + uu(72) - 4*uu(62));
      up(63) - 81 * (uu(62) + uu(64) + uu(53) + uu(73) - 4*uu(63));
      up(64) - 81 * (uu(63) + uu(65) + uu(54) + uu(74) - 4*uu(64));
      up(65) - 81 * (uu(64) + uu(66) + uu(55) + uu(75) - 4*uu(65));
      up(66) - 81 * (uu(65) + uu(67) + uu(56) + uu(76) - 4*uu(66));
      up(67) - 81 * (uu(66) + uu(68) + uu(57) + uu(77) - 4*uu(67));
      up(68) - 81 * (uu(67) + uu(69) + uu(58) + uu(78) - 4*uu(68));
      up(69) - 81 * (uu(68) + uu(60) + uu(59) + uu(79) - 4*uu(69));
      up(72) - 81 * (uu(71) + uu(73) + uu(62) + uu(82) - 4*uu(72));
      up(73) - 81 * (uu(72) + uu(74) + uu(63) + uu(83) - 4*uu(73));
      up(74) - 81 * (uu(73) + uu(75) + uu(64) + uu(84) - 4*uu(74));
      up(75) - 81 * (uu(74) + uu(76) + uu(65) + uu(85) - 4*uu(75));
      up(76) - 81 * (uu(75) + uu(77) + uu(66) + uu(86) - 4*uu(76));
      up(77) - 81 * (uu(76) + uu(78) + uu(67) + uu(87) - 4*uu(77));
      up(78) - 81 * (uu(77) + uu(79) + uu(68) + uu(88) - 4*uu(78));
      up(79) - 81 * (uu(78) + uu(70) + uu(69) + uu(89) - 4*uu(79));
      up(82) - 81 * (uu(81) + uu(83) + uu(72) + uu(92) - 4*uu(82));
      up(83) - 81 * (uu(82) + uu(84) + uu(73) + uu(93) - 4*uu(83));
      up(84) - 81 * (uu(83) + uu(85) + uu(74) + uu(94) - 4*uu(84));
      up(85) - 81 * (uu(84) + uu(86) + uu(75) + uu(95) - 4*uu(85));
      up(86) - 81 * (uu(85) + uu(87) + uu(76) + uu(96) - 4*uu(86));
      up(87) - 81 * (uu(86) + uu(88) + uu(77) + uu(97) - 4*uu(87));
      up(88) - 81 * (uu(87) + uu(89) + uu(78) + uu(98) - 4*uu(88));
      up(89) - 81 * (uu(88) + uu(90) + uu(79) + uu(99) - 4*uu(89));

      tspan = linspace(0, 0.3, 100);
      options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);

      [t, y] = ode15i(@klu, tspan, uu0, up0, options);

      sol = zeros(10, 10, 100);
      for z = 1:100
      for i = 1:10
      sol(:, i, z) = y(z, (((i - 1) * 10) + 1):(i * 10));


      for k = 1:100
      surf(sol(:, :, k))
      axis([0 10 0 10 0 1]);
      title('2-D heat equation semidiscretized to a DAE on the unit square');

      This is the solution at time t = 0.0091:

      Larger problems will be used in order to test the efficiency of the code, because this 2 were solved almost immediately.

      by Francesco Faccio ( at June 20, 2016 01:14 AM

      June 19, 2016

      Amr Mohamed


      Hi all,
      I would like to share my experience through the first weeks of GSoC 2016 with GNU Octave.
      My project is mainly concerned with computational geometry as it aims at creating 2D polygon functions as a part of the geometry package.
      The main function is polybool which permits performing boolean operations on polygons.
      The project’s bitbucket repository can be found here:

      • We started working on the project earlier than the program official start and were able to finish the first required functions for changing the representation of the polygons between cell arrays format and NaN-delimited vectors format.
        These two functions polyjoin/polysplit are implemented as .m files and are used to manipulate the polygons’s representation.
      • Then , we started implementing the ispolycw function that checks the orientation of multiple polygons at the same time using the Boost::Geometry library.
        For self-intersecting polygons, The orientation of the polygon is defined as the orientation of the leftmost point and its two neighbouring points.
      • A Makefile was written to compile and link the written cc files.
        The Makefile was inspired from the sockets package for Octave.
      • Currently, We are focusing on implementing the polybool and we managed to create an initial working version of the function.
        Here are some snapshots for the output of the polybool function:

      theta = linspace(0, 2 * pi, 1000);
      x1 = cos(theta) – 0.5;
      y1 = – sin(theta);
      x2 = x1 + 1;
      y2 = y1 + 1;

      plot(x2, y2)
      [xa, ya] = polybool(‘union’, x1, y1, x2, y2);
      [xb, yb] = polybool(‘intersection’, x2, y2, x1, y1);
      [xc, yc] = polybool(‘xor’, x1, y1, x2, y2);
      [xd, yd] = polybool(‘subtraction’, x1, y1, x2, y2);
      subplot(2, 3, 1)
      axis equal
      patch(xa{1}, ya{1}, ‘FaceColor’, ‘r’)
      axis ([- 2 2 – 2 2])

      subplot(2, 3, 2)
      axis equal
      for k = 1:numel(xb), patch(xb{k}, yb{k}, ‘FaceColor’, ‘r’) end
      axis ([- 2 2 – 2 2])

      subplot(2, 3, 3)
      axis equal
      [x1 y1] = polysplit(x1, y1);
      for k = 1:numel(x1), patch(x1{k}, y1{k}, ‘FaceColor’, ‘b’), end
      axis ([- 2 2 – 2 2])
      title(‘Polygon 1’)

      subplot(2, 3, 4)
      axis equal
      for k = 1:numel(xc), patch(xc{k}, yc{k}, ‘FaceColor’, ‘r’) end
      axis ([- 2 2 – 2 2])

      subplot(2, 3, 5)
      axis equal
      for k = 1:numel(xd), patch(xd{k}, yd{k}, ‘FaceColor’, ‘b’), end
      axis ([- 2 2 – 2 2])

      subplot(2, 3, 6)
      axis equal
      [x2 y2] = polysplit(x2, y2);
      for k = 1:numel(x2), patch(x2{k}, y2{k}, ‘FaceColor’, ‘b’), end
      axis ([- 2 2 – 2 2])
      title(‘Polygon 2’)

      by amrkeleg at June 19, 2016 11:32 PM

      June 17, 2016

      Abhinav Tripathi


      The mid term evaluations of GSoC 2016 are close. This blog post is to summarize the complete work done under the ‘Symbolic Project’ under the organization ‘Octave’.

      I will mention each of the goals that were set for the mid term and then give an overview about what was done for that:

      1a) Octave, Symbolic, and PyTave dev versions installed.
      I built Octave from source successfully on ubuntu 16.04 then used symbolic with it. Also built pytave from source using the dev version of Octave.
      My Symbolic fork on github can be found at:
      My main Pytave fork can be found at:
      I also forked pytave from Colin’s fork (to work on experimental features) which can be found at:

      1b) Some basic communication working between PyTave and Symbolic. This might use some existing implementations in a non-optimal way (e.g., its ok if Symbolic returns objects as xml strings).
      PR #452 added the functionality to convert python types into Octave types including @sym objects using pytave. Then PR #460 fixed the communication when lists/tuples were passed from python to Octave.
      Then we also added proper conversion of tuples and booleans from python to Octave types in pytave repo in the PRs #9 and #11 respectively.

      1c) Most of Symbolic’s tests and doctests continue passing, although some failures ok at this point.
      Tested on Ubuntu 16.04. Most of the tests pass on Linux. There are still many tests which are failing (with the new IPC) but we will work on them once we have a stable pytave IPC mechanism.
      With use of pytave, the errors have now been converted to python exceptions which seems to be the main reason of many tests failing. Also, sometimes it seems that the way we chose to return the outputs is conflicting with how the existing IPC mechanisms were doing it. But, these are minor failures and can be taken care of later as the new IPC matures.

      1d) The above works on at least one OS (probably GNU/Linux).
      It works good on Linux (our local machines).
      We are currently trying to integrate use of pytave on the build-bot. The work can be tracked on the PRs #477 & #478 on github.

      2a) PyTave converts various Python objects into appropriate Octave objects. PyTave needs to be extended to return proper @pyobj when it cannot convert the object. Also, symbolic must be able to convert such objects to @sym objects by calling proper python functions via PyTave (if they are indeed @sym). That is, bypass the current generating of xml strings.
      Instead of editing pytave to convert @sym objects before returning, we incorporated a different mechanism in Symbolic only which converts the objects it gets into @sym objects if necessary using py* functions. The major work was on the PR #452
      Groundwork has been laid for storing objects persistently on python side using @pyobject from pytave. We will also work on improving @pyobject to allow for calling all the attributes of a python object from Octave.

      2b) Improve the BIST test set coverage of both PyTave and Symbolic for any new features added.
      Added some tests to pytave with the PRs #10 and #11 on bitbucket.

      2c) Improve doctest coverage of PyTave.
      Not many doctests were needed so this part is low priority for now. Might still do it this week.

      We also had some Stretch Goals’ which were planned in case we get time and following is the progress for them:

      3) Improve on 2a) with some TBD mechanism: perhaps a pure m-file callback feature in the PyTave export code.
      Since the workaround for @sym conversion was already used in symbolic so this was not required. Moreover, we are working on adding @pyobject support to pytave on Colin’s fork:

      4) The above works on both GNU/Linux and MS Windows.
      Pytave needs to be built on Windows which might take a while. On Linux, we have it working. For building on windows, we are having an extensive and comprehensive discussion and are getting closer to the goal everyday.
      First ‘cygwin’ was tried to build pytave on windows but many tools were not found on cygwin, moreover Octave already has a MSYS environment so we shifted to using something like that. Currently we were trying MSYS2 environment for building pytave, all the libraries were available but we got stuck at a point when trying to link with octave libraries. Octave was build with GCC 4 while MSYS2 has GCC 5 and due to some ABI incompatibility of GCC between the two versions, it was not possible to continue the build.
      So, now we will try to build it from within octave, using the the MSYS bash shell. But, now we need to (probably) build libboost-dev and python-dev using that shell and then move to building pytave with it.

      5) Objects passed from Symbolic back to Python should use PyTave mechanisms.
      The work can be seen in the PR #465 which tries to call py* functions to store the variables into python. We have written a function which stores the required variables to a python list. Currently, it is not very stable and causes octave to crash in some tests but we are working on improving it. It now supports passing lists and 2-D matrices. But, it needs to be extended to support other types like cell-arrays and n-D matrices.

      In the start of the project we had some minor enhancements to Symbolic. In the mid way we also had some enhancements on the pytave side which were needed to move ahead with the aforementioned goals.
      All of my contributions to the repos can be tracked at following links:
      My contributions to octsympy‘s main repo on github
      My contributions to pytave‘s main repo on bitbucket
      My Contribution to Colin’s pytave fork

      This was the summary as per the goals set. Now, a quick summary as per the timeline we decided just to showcase how closely we are following the path we decided at the start of the project:

      30th May — Work on adding the sym conversion to PyTave and cleaning up the conversion mechanism in Symbolic.
      Till 30th May we had already added the mechanism to convert sym objects that come from pytave. We still have to work on it for some cases, as we have many failing tests. But, we have a basic structure that doesn’t rely on any XML stuff for object conversion.

      15th June — Improve tests and doctests. Work on building PyTave and testing on Windows. No more crufty XML!
      Till 15th June we had improved BISTs of pytave. The building of pytave has been underway and we have got rid of the XML stuff in pytave IPC.

      27th June —  Try to get working PyTave with Symbolic on Windows (if needed, use cygwin) 
      We have been following the timeline from the start. We also have other smaller goals in progress aside to this (including build pytave on travis and using pytave to pass variables to python). We have used MSYS2 to build pytave on windows. But that led to some conflict in gcc versions of MSYS2 and the ones present in Octave. Finally it seems that boost and python-dev have to be build using the MSYS environment of Octave and then only pytave can be built using those tools.


      All of the aforementioned goals and timeline can be found at my wiki page:


      For any feedback/comment, feel free to post a comment here…

      by genuinelucifer at June 17, 2016 07:02 PM

      June 15, 2016

      Chiara Segala

      Week 2-3: general exponential schemes

      During these two weeks I implemented the two schemes for a general exponential Runge-Kutta and Rosenbrock integrator, see the schemes described in my second post.
      I also implemented another file for exponential Runge-Kutta integrators based on the following scheme.
      Given the problem
      u'(t) = F(t,u) = A u(t) + g(t, u(t)),
      u(t_0) = u_0,
      the new scheme is
      U_{ni} = u_n + c_i h_n \varphi_1(c_i h_n A) F(t_n,u_n) + h_n \sum_{j=2}^{i-1} a_{ij}(h_n A) D_{nj} ,
      u_{n+1} = u_n + h_n \varphi_1(h_n A) F(t_n,u_n) + h_n \sum_{i=2}^s b_i(h_n A) D_{ni} ,
      D_{nj} = g(t_n + c_j h_n, U_{nj}) - g(t_n, u_n).
      The main motivation for this reformulation is that the vectors $D_{nj}$ are expected to be small in norm so the application of matrix functions to these vectors are more efficient.

      I applied the general pattern to the third-order exponential Runge-Kutta method with tableau

      where $\varphi_{j,k} =\varphi_j (c_k h A) $ and $\varphi_j =\varphi_j (h A) $.

      Then I also applied the general Rosenbrock pattern to the fourth-order exponential method with a third-order error estimator. Its coefficients are

      where $\varphi_j =\varphi_j (h J_n) $.

      I tested the correctness and order of the schemes with the following example, a semilinear parabolic problem
      \frac{\partial U}{\partial t} (x,t) - \frac{\partial^2 U}{\partial x^2} (x,t) = \frac{1}{1 + U(x,t)^2} + \Phi (x,t)
      for $x \in [0,1]$ and $t \in [0,1]$ with homogeneous Dirichlet boundary conditions. $\Phi$ is chosen in such a way that the exact solution of the problem is $U(x,t) = x(1-x)e^{t}$. I discretize this problem in space by standard finite differences with 200 grid points.

      For details see [HO 10] and [HO 05].

      Within this week I will insert on Bitbucket the codes of phi functions and the three general schemes .

      by Chiara Segala ( at June 15, 2016 02:32 AM

      June 14, 2016

      Barbara Lócsi

      Preliminary balancing


      I have created a repository, you can take a look at my work here:

      Preliminary balancing

      Currently in octave eig uses preliminary balancing and it can’t be turned off, while in Matlab[1] it can be.

      §  e = eig(A)
      §  [V,D] = eig(A)
      §  [V,D,W] = eig(A)
      §  e = eig(A,B)
      §  [V,D] = eig(A,B)
      §  [V,D,W] = eig(A,B)
      §  [___] = eig(A,balanceOption)
      §  [___] = eig(A,B,algorithm)
      §  [___] = eig(___,eigvalOption)
      §  lambda = eig (A)
      §  lambda = eig (A, B)
      §  [V, lambda] = eig (A)
      §  [V, lambda] = eig (A, B)

      The aim of this task was to change the *geev LAPACK calls to the extended *geevx which allows us to turn off the balancing. The ability to do this is important becuase the balancing can be harmful:

      “Balancing sometimes seriously degrades accuracy. In particular, one should not balance a matrix after it has been transformed to Hessenberg form. However, we must emphasize that balancing is usually not harmful and often very beneficial. When in doubt, balance.” [2]


      To test the preliminary balancing I used Matlab’s example for „Eigenvalues of a Matrix Whose Elements Differ Dramatically in Scale” [3]

      A = [ 3.0     -2.0      -0.9     2*eps;
           -2.0      4.0       1.0    -eps;
           -eps/4    eps/2    -1.0     0;
           -0.5     -0.5       0.1     1.0];
      [VN,DN] = eig(A,'balance');
      A*VN - VN*DN

      But regardless of the balancing is on or off the results were the same.

        -4.4409e-16   2.2204e-16  -2.9055e-16  -1.6653e-16
         8.8818e-16   2.7756e-16   3.1094e-17  -1.3878e-16
         1.7176e-17   1.5361e-18   6.6297e-18   0.0000e+00
        -6.9389e-17  -4.4409e-16   2.2204e-16   6.9389e-17

      (Which is the more accurate one, the one we want)
      It seems like the matrix is not balanced. Moreover balance(A) is the same as A:


      ans =
         0   0   0   0
         0   0   0   0
         0   0   0   0
         0   0   0   0

      The reason for this is that the stopping criteria for balancing was changed in LAPACK  3.5.0 [4], so it deals more intelligently with cases where the balancing causes loss of accuracy. [5][6][7]

      “However, for the case where A is dense and poorly scaled, the new algorithm will still balance the matrix and improve the eigenvalue condition number. If accurate eigenvectors are desired, then one should consider not balancing the matrix.” [5]

       Other matrix

      A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1]; % drop the eps terms for now

      Balancing this matrix is not harmful, so LAPACK will balance as we expect it to do.
      A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1]; % drop the eps terms for now
      [VN,DN] = eig(A,'nobalance');
      A*VN - VN*DN
      [VN,DN] = eig(A,'balance');
      A*VN - VN*DN
      [VN,DN] = eig(A);
      A*VN - VN*DN

      ans =

        -4.4409e-16   2.2204e-16  -2.4947e-16  -1.6653e-16
         8.8818e-16   3.8858e-16  -4.1531e-17   1.3878e-16
         0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
        -8.3267e-17  -4.4409e-16   2.2204e-16  -4.1633e-17

      ans =

         0.0000e+00  -4.4409e-16   2.2204e-16   1.6653e-16
         0.0000e+00   0.0000e+00  -2.2204e-16  -1.3878e-16
         0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
         0.0000e+00   1.3878e-17   0.0000e+00   0.0000e+00

      ans =

         0.0000e+00  -4.4409e-16   2.2204e-16   1.6653e-16
         0.0000e+00   0.0000e+00  -2.2204e-16  -1.3878e-16
         0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
         0.0000e+00   1.3878e-17   0.0000e+00   0.0000e+00

      On this example we can see that the eig(A) is working as the eig(A,'balance') which is a behaviour we want.

      So the ability to turn off the balancing is not as important as it was before LAPACK 3.5.0. but there are some cases when the balancing is bad, so the ability to turn off the balancing is still important.

      “However, for the case where A is dense and poorly scaled, the new algorithm will still balance the matrix and improve the eigenvalue condition number. If accurate eigenvectors are desired, then one should consider not balancing the matrix.” [5]

      by Barbara Lócsi ( at June 14, 2016 05:18 AM

      June 09, 2016

      Francesco Faccio

      First Goals and next steps

      MathJax TeX Test PageHello!

      I have not written a post for a while because I have had some health issues.

      During the first two weeks of GSoC I have worked on Autotools and I have compiled Octave with link to SUNDIALS. The first step for doing this was to check the presence and usability of ida.h in, so I used the macro OCTAVE_CHECK_LIB which also sets the flags CPPFLAGS, LDFLAGS and LIBS. Then I set the right configuration variables in the build-aux folder and modified build-env namespace. Finally I wrote a dld-function which includes ida.h and calls the function IDACreate from SUNDIALS which returns a pointer to the IDA memory structure.
      This dld-function generates an oct-file which can be executed from Octave.

      All these changes are visible in my public repository on Bitbucket:

      In the next few days I will further investigate the recursive dependencies of SUNDIALS and their license and set up the correct build flags for such dependencies, I will write more tests in in order to check the availability of functions and headers of the library.

      After discussing with mentors we decided to start the implementation of ode15i because it's more close to IDA and more general than ode15s. Once ode15i will be written, ode15s will be built around it.

      We have also decided which are the next steps before midterm evaluation:
      • implement a minimal .oct wrapper for IDA in Octave with a primitive interface such as $[t , y] = ode15i (odefun, tspan, y0, yp0, Jacobian)$
        that invokes IDA with all options set to default values

      • use two benchmark problems to test the correctness and speed of the code:
        I will compare it with the C implementation of SUNDIALS and with the m-file implementation relying on the mex interface of SUNDIALS 

      As benchmark problems we have chosen two examples which deal with dense and sparse methods.

      The first one regards Robertson chemical kinetics problem, in which differential equations are given for species $y_{1}$ and $y_{2}$ while an algebraic equation determines $y_{3}$. The equations for the species concentrations $y_{i}(t)$ are:

      \begin{eqnarray*} \begin{cases} y_{1}^{'} = -0.04y_{1} + 10^{4}y_{2}y_{3} \\ y_{2}^{'} = 0.04y_{1} - 10^{4}y_{2}y_{3} - 3\cdot 10^{7}y_{2}^{2} \\ 0 = y_{1} + y_{2} + y_{3} - 1 \end{cases} \end{eqnarray*}

      The initial values are taken as $y_{1} = 1$, $y_{2} = 0$ and $y_{3} = 0$. This example computes the three concentration components on the interval from $t = 0$ through $t = 4\cdot 10^{10}$.

      This is the plot of the solution (the value of $y_{2}$ is multiplied by a factor of $10^{4}$).

      Dense methods of IDA are applied for solving this problem.

      The second problem is a $2D$ heat equation, semidiscretized to a DAE. The DAE system arises from the Dirichlet boundary condition $u = 0$, along with the differential equations arising from the discretization of the interior of the region.
      The domain is the unit square $\Omega = \{0 \leq x, y \geq 1\}$ and the equations solved are:

      \begin{eqnarray*} \begin{cases} \partial u/\partial t = u_{xx} + u_{yy} & (x, y) \in \Omega \\ u = 0 & (x, y) \in \partial \Omega \end{cases} \end{eqnarray*}

      The time interval is $0 \leq t \leq 10.24$, and the initial conditions are $u = 16x(1 − x)y(1 − y)$.
      We discretize the PDE system (plus boundary conditions) with central differencing on a $10 \times 10$ mesh, so as to obtain a DAE system of size $N = 100$. The dependent variable vector $u$ consists of the values $u(x_{j}, y_{k}, t)$ grouped first by $x$, and then by $y$. Each discrete boundary condition becomes an algebraic equation within the DAE system.

      In this problem IDA's sparse direct methods are used and the Jacobian is stored in compressed sparse column (CSC) format.

      Regarding functions which deal with input ode check:
      Functions check_input and set_ode_options, which I started to write before the beginning of the coding period, will be improved after the midterm evaluation.

      by Francesco Faccio ( at June 09, 2016 02:49 PM

      June 08, 2016

      Cristiano Dorigo

      Presentation of the Project and Timeline

      Dear all,

      I'm Cristiano Dorigo, a bachelor graduated student in Applied Mathematics in Verona, Italy, and now I'm a student of Master Degree in Mathematics in the University of Verona.

      I was selected as partecipant student for the SOCIS 2016 under GNU Octave for the project "Improve iterative methods for sparse linear systems".

      The idea of this project is at first to improve the methods actually implemented in Octave to solve linear systems as gmres, pcg, bicg, bicgstab, qmr, pcr. Ideed these algorithms have, for example, a similar structure in the initial part (where the inputs are checked). Other possible improvement can be the synchronization of their error messages or to check the documentation. Another critical point is to check the compatibility of the results with Matlab, i.e. to check, for example, if the quality of the solution with the same method is equal both in Matlab and Octave, to see if they use the same number of iterations, same residual,...

      The second step is to implement, as much as possible, some iterative methods that, for example, are in Matlab, but not in Octave. The suggested methods are minres, symmlq, tfqmr, bicgstabl, lsqr. The order of the implementation can be decided with the mentors giving priority to the useful or the necessary ones.

      I think to split the work in the same two step mentioned before, i.e. improvements and implementation.

      • The first two weeks I want to study the methods already implemented but that I haven't faced off during my student career, i.e. bicg, bicgstab, qmr and pcr and to study their codes.
      • The second two/three weeks I want at first to consult the mentors to decide where and how to improve these methods, then to apply these ideas.
      • The remaining time will be spent to implement the new methods. I want to implement these methods one by one, after a discussion with the mentors about which methods have priority.

      The reference book is "Iterative Methods for sparse linear systems" 2-nd Edition by Yousef Saad, which is available here

      by Cristiano Dorigo ( at June 08, 2016 03:08 AM

      June 01, 2016

      Chiara Segala

      Week 1: phi functions

      pre.cjk { font-family: "Nimbus Mono L",monospace; }p { margin-bottom: 0.25cm; line-height: 120%; }

      During the first week of GSoC, I wrote four m-files for matrix functions, phi1m,..., phi4m.
      I implemented the fourfunctions, based on [BSW 07].
      Matrix functions are defined by the recurrence relation

      In my files, I use an algorithm based on first a variant of the scaling and squaring approach and after scaling, I use a Padé approximation.
      Below the code phi1m.m

      function [N, PHI0] = phi1m (A, d)

      Nis φ1(A) and PHIOisφ0(A). dis the order of the Padé approximation.
      First, I scale A by power of 2 so that its norm is < ½


      s= min (ceil (max (0,1+log2 (norm (A, inf)))), 1023);

      A = A/2^s;

      then I use a (d,d)-Padé approximation

      where the polynomials are

      I write polynomials in Horner form.

      ID = eye (size (A));


      Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

      Dcoeff = ((-1)^(i))*prod ([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

      N= Ncoeff;

      D= Dcoeff;

      for i = (d-1):-1:0

         Ncoeff = sum (cumprod ([1, d-(0:i-1)])./cumprod ([1, 2*d+l-(0:i-1)]).*(-1).^(0:i)./(factorial ((0:i)).*factorial (l+i-(0:i))));

         N = A*N + Ncoeff * ID ;

         Dcoeff = ((-1)^(i))* prod([(d-i+1):1:d])/(prod ([(2*d+l-i+1):1:(2*d+l)])*factorial (i));

         D = A*D + Dcoeff * ID ;


      N = full (D\N);

      and finally I use the scaling relations

      PHI0 = A*N+ID;

      for i = 1:s

         N = (PHI0+ID)*N/2;

         PHI0 = PHI0*PHI0;


      in the same way I wrote the other three files.  As soon as possible, I will create a repository on Bitbucket and I will put the codes there.

      by Chiara Segala ( at June 01, 2016 04:25 AM

      Phi functions and general exponential schemes

      p { margin-bottom: 0.25cm; line-height: 120%; }

      May 23 will start the coding period and I'm trying to figure out how I will implement my functions before the mid-term evaluation.

      With the advice of my mentors Marco and Jacopo, I decided to start with the implementation of the phi functions, necessary to calculate the matrix functions in the two schemes for a general exponential Runge-Kutta and Rosenbrock integrator.
      These schemes will not be really fast and efficient, but I will use them as a reference when I go to implement the official methods. It will be useful to verify the correctness of my codes.

      As regards the implementation of the phi functions I will refer to

      [BSW 07] “EXPINT — A MATLAB Package for Exponential Integrators”, Havard Berland, Bard Skaflestad and Will M. Wright, 2007,
      DOI: 10.1145/1206040.1206044, webpage (software without a license).

      While the general schemes that then I'm going to implement are as follows:

      • Exponential Runge-Kutta integrators
      p { margin-bottom: 0.25cm; line-height: 120%; }
      Consider a problem of the form



      p { margin-bottom: 0.25cm; line-height: 120%; }

      the numerical exponential Runge Kutta scheme for its solution is


      and the coefficients and are constructed from exponential functions or approximations of such functions.

      • Exponential Rosenbrock integrators
      Consider a problem of the form


      the numerical exponential Rosenbrockscheme for its solution is


      for details about formulas see [HO 10].

      by Chiara Segala ( at June 01, 2016 04:14 AM

      May 31, 2016

      Abhinav Tripathi


      So, according to the timeline that we decided, the first phase has been complete. The timeline stated:

      30th May — Work on adding the sym conversion to PyTave and cleaning up the conversion mechanism in Symbolic.

      I think we have done some good work to fulfill the objective. We have had many changes to add proper @sym conversion on Symbolic side. The following 2 PRs take care of that:

      Also, we worked on the PyTave side to fix conversion of tuples and booleans properly from python to octave. It can bee seen in the following PRs:


      As far as the midterm goals are concerned. We have completed 1a), 1b) 1d) & 2a).

      Most of 1c) is also done. For now it is acceptable is my view. We now move to 2b) and 2c).

      The next goal (as per timeline is)-
      15th June — Improve tests and doctests. Work on building PyTave and testing on Windows. No more crufty XML!

      We already got rid of the XML stuff in the new PyTave IPC. We will now be working on adding and improving test coverage of  PyTave (both BISTs and doctests). In approximately 10 days,  we will move on to start trying to build PyTave on windows.

      by genuinelucifer at May 31, 2016 06:27 PM


      We (me and both mentors) decided upon the following timeline to be followed during the course of this project:

      30th May — Work on adding the sym conversion to PyTave and cleaning up the conversion mechanism in Symbolic.
      15th June — Improve tests and doctests. Work on building PyTave and testing on Windows. No more crufty XML!
      27th June — Try to get working PyTave with Symbolic on Windows (if needed, use cygwin) [Mid Term Evaluations]

      5th July — Get a successfully working Symbolic with PyTave (all the tests and features working). Continue work on Goal 3.
      20th July — Finalize implementation for Goal 3.
      31st July — Work on improvements to PyTave such as adding support for Python objects from within Octave (printing, calling methods, saving etc…).
      10th August — Recode @sym methods as required to take benefit of PyTave. Also, add some other methods into Symbolic from #215
      15th August — “Blue-skies” stuff. Try to fix the unicode utf8 rendering on Windows. Explore possibility of incorporating PyTave into Octave. [gsoc final code submission begins]
      23rd August — Finish all the code, test on buildbots and Windows. Submit to Google. [Final deadline for code submission]
      Afterwards — Complete any goals which are left. Continue contributing to Octave…


      Following goals were set for the mid term evaluations :

      1a) Octave, Symbolic, and PyTave dev versions installed.

      1b) Some basic communication working between PyTave and Symbolic. This might use some existing implementations in a non-optimal way (e.g., its ok if Symbolic returns objects as xml strings).

      1c) Most of Symbolic’s tests and doctests continue passing, although some failures ok at this point.

      1d) The above works on at least one OS (probably GNU/Linux).

      2a) PyTave converts various Python objects into appropriate Octave objects. PyTave needs to be extended to return proper @pyobj when it cannot convert the object. Also, symbolic must be able to convert such objects to @sym objects by calling proper python functions via PyTave (if they are indeed @sym). That is, bypass the current generating of xml strings.

      2b) Improve the BIST test set coverage of both PyTave and Symbolic for any new features added.

      2c) Improve doctest coverage of PyTave.

      Stretch Goals: 

      3) Improve on 2a) with some TBD mechanism: perhaps a pure m-file callback feature in the PyTave export code.

      4) The above works on both GNU/Linux and MS Windows.

      5) Objects passed from Symbolic back to Python should use PyTave mechanisms.

      by genuinelucifer at May 31, 2016 05:59 PM

      May 19, 2016

      Amr Mohamed


      It’s only four days for GSoC coding period to start and the project seems to be in the right track.
      During the last few weeks, I have finished the first .m files for the polyplit and polyjoin functions.
      Currently i am trying to prepare a makefile for compiling the rest of the functions as they will be .cc files not .m files.

      by amrkeleg at May 19, 2016 09:00 PM

      May 17, 2016

      Chiara Segala

      Presentation and timeline


      I'm Chiara Segala, graduated at University of Verona, Italy. I attended my bachelor's degree in applied mathematics and I am doing now the second year of the master's degree.
      I was selected for the GSoC 2016 with the project Exponential Integrators (GNU Octave organization).
      Exponential integrators are a class of numerical methods for the solution of partial and ordinary
      differential equations.

      This is an estimated TIMELINE of my work:

      05/05 - 23/05
      Study the theory about the exponential Runge-Kutta and Rosenbrock-type integrators:
      [HO 10] “Exponential integrators”, Marlis Hochbruck and Alexander Ostermann, 2010, DOI: 10.1017/S0962492910000048
      [HO 05] "Explicit exponential Runge-Kutta methods for semilinear parabolic problems", Marlis Hochbruck and Alexander Ostermann, 2005, DOI: 10.1137/040611434, preprint.
      Familiarize with the other ODE solvers in Octave, odepkg.
      Take a look at
      [J 14] “EXPODE - Advanced Exponential Time Integration Toolbox for MATLAB”, Georg Jansing, 2014, webpage.
      Study the theory and the features of the expmv code:
      [HAM 11] “Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators”, Awad H. Al-Mohy and Nicholas J. Higham, 2011, DOI: 10.1137/100788860.

      23/05 - 27/06
      week 1 : implementation of phi functions
      week 2 : implementation of a scheme for a general exponential Runge-Kutta integrator, see section 2.3 of [HO 10], using the phi functions
      week 3 : implementation of a scheme for a general exponential Rosenbrock integrator , see section 2.4 of [HO 10], using the phi functions
      week 4 : implementation of a method for the construction of matrix Ã, see theorem 2.1 of [HAM 11]

      Mid-term evaluation

      27/06 - 15/08
      week 1-2 : implementation of an advanced method with adaptive time stepping from Runge-Kutta family, using expmv
      week 3-4 : implementation of an advanced method with adaptive time stepping from Rosenbrock-type family, using expmv
      week 5-6 : implementation of validation tests, e.g. analyze the order of convergence of the methods and some examples
      week 7 : add some improvement and code clean up
      week 8 : write documentation

      15/08 - 22/08 : review of all the work

      by Chiara Segala ( at May 17, 2016 05:39 AM

      May 15, 2016

      Francesco Faccio



      This is a Timeline for the project ode15s.

      As discussed with the mentors, our goals for the mid-term evaluations are to build Octave with all the dependencies of SUNDIALS and to create an m-file which deals with the input parameters and the options of a generic ODE/DAE solver.
      The final goal is  to have a well tested and documented implementation of ode15s.

      Community Bonding Period: 
      -Familiarize with Autotools and the structure of Octave
      -Study the documentation of SUNDIALS and Oct-files

      Week 1-2 (May 23 - Jun 5):
      -Add SUNDIALS as a dependency and build Octave from source (I will usa a dld_function which calls a function of SUNDIALS)

      Week 2-3 (Jun 6 - 19):
      -Write an m-file which deals with the input of a generic ode/dae solver

      Midterm Evaluations

      Week 4 (Jun 20 - 26):
      -Design the code of ode15s (I will choose which functions of SUNDIALS will be used)

      Week 5-6 (Jun 27 - Jul 10):
      -Write Oct-files

      Week 7 (Jul 11 - 17): 
      -Write tests 

      Week 8-9 (Jul 18 - 31):
      -Test compatibility between Matlab and Octave
      -Test the performance of the algorithm

      Week 10-11 (Aug 1 - 14):
      -Write the documentation and perform more tests

      Week 12 (Aug 15 - 21)
      -Review of the work

      OPTIONAL: If I finish the work early, I will try to write a (slower) version of ode15s which uses DASPK or DASSL (this is for people who don't have SUNDIALS installed)

      Final Evaluations

      by Francesco Faccio ( at May 15, 2016 02:29 PM