domenica 24 novembre 2013

Backward Euler with Inexact Newton solver

Backward Euler

The backward Euler method is one of the most common methods used to solve ODEs. It is an implicit method of order one but its strenght lies in its A-stability property.

Given a system of ODEs of the form: $$\frac{d\mathbf{y}}{dt} = \mathbf{f}(t,\mathbf{y}) \\ \mathbf{y}(t_0) = \mathbf{y}_0\ ,$$ the backward Euler method is written as: $$ \mathbf{y}_{k+1} = \mathbf{y}_k + \Delta t\ \mathbf{f}(t_{k+1},\mathbf{y}_{k+1})\ .$$

If we define u to be the exact solution then we have that: $$ \frac{d\mathbf{u}}{dt} = \mathbf{f}(t,\mathbf{u}) \\ \mathbf{u}(t_k) = \mathbf{y}_k $$ So that $$\mathbf u(t_{k+1}) = \mathbf u(t_k) + \int_{t_k}^{t_{k+1}} \mathbf{f}(t,\mathbf u) dt\ ,$$ which, under some conditions of regularity, can be rewritten as $$\mathbf u(t_{k+1}) = \Delta t\ \mathbf f(t_{k+1},\mathbf u(t_{k+1})) - \frac{{\Delta t}^2}{2} \mathbf f'(\xi)\ .$$

From these results it can be derived that the local truncation error (LTE) is $$ \mathbf e_{k+1} = \mathbf y_{k+1} - \mathbf u(t_{k+1}) = \mathbf o({\Delta t}^2)\ .$$

Inexact Newton

Now we have to solve the non-linear system built up with the backward Euler method. We want to use an inexact Newton solver. The classical Newton method is written as: given a non-linear functions system F(x) = 0 with a given initial condition, solve the linearized system: $$ \mathbf F'(\mathbf x_k)\mathbf s_k + \mathbf F(\mathbf x_k) = \mathbf 0 \\ \mathbf x_{k+1} = \mathbf x_k + \mathbf s_k\ .$$ Go on iterating this method until a given tolerance is reached.

In many cases, especially when we don't have an explicit expression for the Jacobian or it can't be inverted, we can use an inexact Newton method: $$ || \mathbf F'(\mathbf x_k)\mathbf s_k + \mathbf F(\mathbf x_k) || \leq \eta_k || \mathbf F(\mathbf x_k)||\ ,$$ where \(\eta_k\) is said to be the forcing term. The linearized system can be solved with an iterative method such as GMRES or preconditioned conjugate gradient.

There are two possible choices for the forcing term:

  1. first choice $$ \eta_k = \frac{|| \mathbf F(\mathbf x_k)-\mathbf F(\mathbf x_{k-1}) -\mathbf F'(\mathbf x_{k-1})\mathbf s_{k-1} ||}{||\mathbf F(\mathbf x_{k-1}) ||}\ ;$$
  2. second choice $$ \eta_k = \gamma \Big(\frac{|| \mathbf F(\mathbf x_k) ||}{|| \mathbf F(\mathbf x_{k-1}) ||}\Big)^{\alpha}\ .$$

Choosing the forcing term with one of the two previous possibilities, rather than imposing a fixed tolerance, it is possible to avoid the phenomenon of oversolving. Another advantage is that if we use the GMRES as the linear solver it's not necessary to know the Jacobian nor to datermine it, because we can approximate the term $$ \mathbf F'(\mathbf x_{k})\mathbf s_k $$ with a first order finite differences formulae: $$ \mathbf F'(\mathbf x_k)\mathbf s_k \approx \frac{\mathbf F(\mathbf x_k + \delta\mathbf s_k) - \mathbf F(\mathbf x_k)}{\delta} \ .$$

The signature of inexact_newton is the same of fsolve:

[x,fval,exitflag,output,jacobian]=inexact_newton(fun,x0,options)

venerdì 15 novembre 2013

Ode solvers options

Ode solvers options

What I want to do at this moment of the project is to add to all the solvers I've written since now, the ability to manage all the possible options of the ODE structure (defined previously in odeset structure) that are already managed in corresponding MATLAB ode solvers; so that the solvers I've written will be totally MATLAB-compatible and this will be the same also for the new solvers I'm going to implement.

mercoledì 6 novembre 2013

Timeline

Timeline

The following list explains the main points of the timeline for my SOCIS project:

  1. adapt odeset and odeget from Octave odepkg to be completely MATLAB-compatible: all option names supported by MATLAB should be available; furthermore we want the possibility to add new option names that will be useful for example for symplectic integrators;
  2. implement a new code structure for one-step integrators ode12, ode23, ode45, odefwe (forward Euler), odebwe (backward Euler) and inexact solvers, subdividing the execution into 'blocks' in order to be able to do an easier optimization and an easier extension to new functionalities; in particular implement the steppers which will be the functions called to do one integration step and then implement the integrate functions which will integrate, in different ways, over the whole time interval, calling the respective steppers. Afterwards adapt the interface of ode45 to be completely MATLAB-compatible;
  3. implement MATLAB-compatible version of deval;
  4. implement the following geometric integrators as new steppers:
    • symplectic Euler;
    • Stormer-Verlet;
    • velocity-Verlet;
    • spectral variational integrators;
    • SHAKE;
    • RATTLE.
  5. optimization, error checking and documentation of all the code written.


At the moment point 1 is completed; in point 2 I still have to adapt the interface of ode45 and to refine the implemented code. Point 3 has to be done. In point 4 I have implemented sympelctic Euler, Stormer-Verlet and spectral variational integrators but these need to be adapted to the code structure 'steppers', 'integrate functions'; the other integrators must be implemented. Finally, during the implementation I will start to write a bit of documentation and at the end I will do optimizations and debugging.

sabato 2 novembre 2013

Odeset & odeget

Odeset & odeget

Odeset and odeget functions allow to build a valid ODE options structure. They are already available in Octave odepkg but they are not perfectly compatible with MATLAB odeset and odeget functions. Furthermore, for geometric integrators like Spectral Variational Integrators,
I will need new options into the ODE options structure, which now are not admitted in Octave odepkg.

So that I have written my own versions of odeget.m, odeset.m and the new function ode_struct_value_check.m in order to have full compatibility with MATLAB odeset and odeget (their same behaviour on the basis of their official documentation) and also the possibility to add new field names which i will need for future solvers of this project.

The new fields are the union of MATLAB ODE options and Octave ODE options, plus my new options (default values are in square brackets):

'             AbsTol:  scalar or vector, >0, [1.e-6]'
'                BDF:  binary, {on, [off]}'
'             Events:  '
'           Explicit:  binary, {yes, no, []}'
'      InexactSolver:  switch, {inexact_newton, []}'
'       InitialSlope:  vector, []'
'        InitialStep:  scalar, >0, []'
'           Jacobian:  matrix or function_handle, []'
'          JConstant:  binary, {on, [off]}'
'           JPattern:  sparse matrix, []'
'               Mass:  matrix or function_handle, []'
'       MassConstant:  binary, {on, [off]}'
'       MassSingular:  switch, {yes, [maybe], no}'
'MaxNewtonIterations:  scalar, integer, >0, [1.e3]'
'           MaxOrder:  switch, {1, 2, 3, 4, [5]}'
'            MaxStep:  scalar, >0, []'
'   MStateDependence:  switch, {none, [weak], strong}'
'          MvPattern:  sparse matrix, []'
'          NewtonTol:  scalar or vector, >0, []'
'        NonNegative:  vector of integers, []'
'        NormControl:  binary, {on, [off]}'
'          OutputFcn:  function_handle, []'
'         OutputSave:  scalar, integer, >0, []'
'          OutputSel:  scalar or vector, []'
'   PolynomialDegree:  scalar or vector, integer, >0, []'
'    QuadratureOrder:  scalar or vector, integer, >0, []'
'             Refine:  scalar, integer, >0, []'
'             RelTol:  scalar, >0, [1.e-3]'
'              Stats:  binary, {on, [off]}'
'     TimeStepNumber:  scalar, integer, >0, []'
'       TimeStepSize:  scalar, >0, []'
'         Vectorized:  binary, {on, [off]}'

The usage of odeset will be one of the following:

odeset
opt = odeset()
opt = odeset(old_ODE_options,new_ODE_options)
opt = odeset(old_ODE_options,'field1',value1,'field2',value2,...)
opt = odeset('field1',value1,'field2',value2,'field3',value3,...)

The usage of odeget will be one of the following:

option_value = odeget(ODE_structure,'field')
option_value = odeget(ODE_structure,'field',default_value)
option_value = odeget(any_struct,'field',default_value,'fast')
The last usage is needed for MATLAB compatibility and represents a fast access to the given structure with no error checking on options names.

Fuzzy string comparison

Users may do mistakes while typing ODE option names or they may want to write everything in low cases for a faster coding. So that my odeset and odeget functions behave in a user-friendly way, calling the fuzzy_compare function which determines the distance of each structure option name from the given word and returns the indices of those options whose distance is under a given tolerance.

The idea is that word cases are not relevant, if we want an exact matching then the tolerance will be 0, if no index is returned then it will definitely be a typing error, if one index is returned but words don't match exactly a warning will be displayed saying that the program is going on assuming that the user was intending the closest option name, otherwise all the fields whose distance is under the tolerance will be displayed and the user will be asked to insert the name again.

Signature of this function follows.

res = fuzzy_compare(string1,string_set,correctness)
  • string1 must be a string and is the option name we're looking for;
  • string_set must be a cell array of strings or a column vector of strings; it represents the set of option names;
  • correctness is the tolerance we want. It is an optional input argument;
We may set correctness equal to 'exact' or 0 and in this case fuzzy_compare will look only for exact matching; if we set correctness to any positive integer it will be the tolerance of the method (that is the maximum distance, from the given string, accepted). If we don't specify anything for correctness then the tolerance will be calculated as a function of the minimum distance and of the length of the given string: the less the distance, the less the tolerance; the greater the length, the greater the tolerance.
tolerance = minimus + floor(length(string1)/4)*floor(minimus/3)
There exist many definitions of distance between strings. I've chosen the Levenshtein distance.

Levenshtein distance

The Levenshtein distance is a string metric and is equal to the minimum number of single-characters edit (insertion, deletion and substitution) required to change one word into the other. This is the main algorithm of levenshtein function:

function [dist,d] = levenshtein(string1,string2)
  m = length(string1);
  n = length(string2);

  d = zeros(m+1,n+1);
  d(2:m+1,1) = [1:1:m];
  d(1,2:n+1) = [1:1:n];

  for j=2:1:n+1
    for k=2:1:m+1
      if(string1(k-1)==string2(j-1))
        d(k,j)=d(k-1,j-1);
      else
        d(k,j)=min(d(k-1,j)+1,min(d(k,j-1)+1,d(k-1,j-1)+1));
      end
    end
  end

  dist = d(m+1,n+1);
end

All the code is available under the terms of GNU General Public License at my_octave-odepkg.