## Copyright (C) 2013 Simone Panza
## Copyright (C) 2013 Fabio Cisaria
## Copyright (C) 2013 Carlo de Falco
##
## This file is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published
## by the Free Software Foundation; either version 3 of the License,
## or (at your option) any later version.
##
## This file is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this file; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {[@var{t}, @var{y}] =} odebwenk @
## (@var{f}, @var{tspan}, @var{y0}, @var{options}, @var{solver}) 
## Solve a system of ODEs by means of Backward Euler method using the
## inexact Newton method to solve nonlinear  
## iterations and a Krylov method to solve linear algebraic systems. 
## Step adaptation is implemented estimating the local truncation
## error by comparison with a second order predictor-corrector 
## method.
##
## Input:
## @itemize @minus
## @item @var{f} is a function handle for the system of ODEs, in the
## normal form dy/dt = f(t,y)  
## @item @var{tspan} is a 2 element vector containing initial and
## final time  
## @item @var{y0} is the initial condition y(t=0)
## @item @var{options} option structure creted by odeset. Supported
## options are: 
## RelTol, AbsTol, NewtonTol, Mass 
## @item @var{solver} is the name of the iterative linear solver
## (default GMRES) 
## @item @var{restart} in case the selected solver is GMRES, sets the
## restart parameter. otherwise it is uneffective 
## @end itemize 
##
## Output:
## @itemize @minus
## @item @var{tout} is th vector of times at which the solution is given
## @item @var{yout} is a matrix such that @code{size (yout) =
## [numel(y0), numel(tout)]}  
## @end itemize
## @seealso{gmres, pcg, bicgstab, cgs, inexact_newton}
## @end deftypefn

## Authors Simone Panza and Fabio Cisaria
## v1.1
## date 11-03-2013

function [tout, yout] = odebwenk
         odebwenk  (f, tspan, y0, options, solver = [], restart = 20)

  ## Initialization
  if (nargin < 4)
    options = odeset (); 
  endif
  
  ## user-defined relative tolerance
  if (isfield (options, 'RelTol') && ! isempty (options.RelTol))
    rtol = options.RelTol;
  else
    rtol = 1.0e-3;
  endif

  ## user-defined absolute tolerance
  if (isfield (options, 'AbsTol') && ! isempty (options.AbsTol)) 
    atol = options.AbsTol;
  else
    atol = 1.0e-6;
  endif

  ## user-defined Newton abs tolerance
  if (isfield (options, 'NewtonTol') && ! isempty (options.NewtonTol)) 
    newt_tol = options.NewtonTol;
  else
    newt_tol = 1.0e-6;
  endif
  
  ## mass matrix
  if (isfield (options, 'Mass') && ! isempty (options.Mass)) 
    ## mass matrix
    if(isnumeric(options.Mass))
      M = options.Mass; ## mass matrix
      isfunction = false;
    elseif(isa (options.Mass, 'function_handle')) 
      ## if mass matrix is a function_handle 
      isfunction = true;
    endif
  else ## default value for mass matrix : identity matrix
    isfunction = false;
    M = eye (numel (y0));
  endif
  
  t = tspan(1);
  tfinal = tspan(2);

  ## user-defined max step size
  if (isfield (options, 'MaxStep') && ! isempty (options.MaxStep)) 
    hmax = options.MaxStep;
  else
    hmax = .1 * abs (tfinal - t);
  endif

  ## user-defined initial step size
  if (isfield (options, 'InitialStep')
      && ! isempty (options.InitialStep)) 
    h = options.InitialStep;
  else
    h = (tfinal - t) * .05; ## initial guess at a step size
  end 

  hmin = min (16 * eps (tfinal - t), h);
  
  yout = y = y0(:); ## this always creates a column vector, y
  tout = t;         ## first output time

  
  ## arbitrary number of nonlinear iterations per step 
  niter = 100; 
  
  ## MAIN LOOP
  while (t < tfinal) && (h >= hmin)
    
    if t + h > tfinal
      h = tfinal - t; 
    endif
    
    t_new = t + h; ## try a step size h
    
    ## define implicit (nonlinear) function to be solved
    if (isfunction)   ## in the case Mass matrix is dynamic:
      
      ## first order: BE. M = M(t, y) 
      be = @(y_new) M (t_new, y_new) * y + ...
           h * f (t_new, y_new) - M (t_new, y_new) * y_new; 

    else ## constant mass matrix

      ## first order: BE. M = nxn cost
      be = @(y_new) M * y + h * f(t_new, y_new) - M * y_new; 

    endif

    ## solve implicit function by means of a Inexact Newton method
    [y_new_be, yiter] = inexact_newton (be, y, newt_tol, niter,  ...
                                        [], [], solver, restart);
    
    
    ## Second order approximation
    if (isfunction)   ## in the case Mass matrix is dynamic:
      ## CN approximated. M = M(t,y)
      y_new_cn = y + h/2 * (M(t_new, y_new_be) \ ...
                             (f(t,y) + f(t_new, y_new_be))); 
    else ## constant mass matrix
      ## CN approximated. M = nxn cost
      y_new_cn = y + h/2 * (M \ (f(t,y) + f(t_new, y_new_be))); 
    endif
    
    ## o(H^3) estimate for the LTE
    err_est = abs (y_new_be - y_new_cn);
    
    tau = max (rtol .* abs (y), atol);
                                                                                                                                     
    if (all (err_est <= tau)) ## the error lies within the tolerance
      t = t + h; ## take a step			
      y = y_new_be; ## take BE solution
      
      tout = [tout t];	## update time vector
      yout = [yout y];	## update solution vector

      ## adaptively update the step size
      factor = 1.0;

      if (any (err_est))
        
        ## see [1]; .9 as a safety margin;
        ## 1/2 exponent because i am using a 1st order formula 
        factor = .9 * sqrt (min (tau(err_est > 0) ./  ...
                                 err_est(err_est > 0))); 
        
      endif

      h = min (hmax, h * factor); ## take the smallest of the
                                  ## increments: the most conservative

      else ## reject the step, try with a smaller one
      if (h <= hmin)
    	error (["odebwenk: requested step-size too small at t = ", ...
                "%g, h = %g, err = %g"], t, h, norm (err_est, inf)) 
      endif
      
      ## same technique to reduce the step
      ## see [1]; .9 as a safety margin; 1/2 exponent because we're
      ## using a 1st order formula 
      factor = .9 * sqrt (min ((tau(err_est > 0) ./  ...
                                err_est(err_est > 0))(:)));  
      h = max (hmin, h * factor);

    endif
    
    
  endwhile

  if (t < tfinal)
    error (["odebwenk: requested step-size too small at t = %g, h ",
            "= %g, err = %g"], t, h, err) 
  endif
  
endfunction

%!demo
%! tspan = [0 10];
%! opt = odeset ('AbsTol', 1e-6, 'RelTol', 1e-3, 'NewtonTol', 1e-6,
%! 'InitialStep', .1); 
%! N = 600;
%! b = ones (N, 1);
%! A = spdiags ([-b, 2*b, -b]*(N^2), -1:1, N, N);
%! ystat = zeros (N, 1);
%! ystat(2:end-1) = A(2:end-1,2:end-1) \ b(2:end-1);
%! f = @(t, y) [0; (b-(A*y+exp(y)/100))(2:N-1); 0];
%! [t, y] = odebwenk (f, tspan, zeros (N, 1), opt, 'pcg');
%! plot (1:N, ystat, 1:N, y(:, end), 'x')
%! legend ('stationary', 'transient')
