## Copyright (C) 2010 Olaf Till <olaf.till@uni-jena.de>
##
## This program 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 program 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 program; If not, see <http://www.gnu.org/licenses/>.

function ret = odetest ()

  ret.p1.helpers.g = @ p1_g;

  ret.p1.helpers.u = @ p1_u;

  ret.p1.ode = @ (x, t) p1_ode (x, t, ret);

  ret.p1.x0 = @ p1_x0;

  ret.p1.jac = @ (x, t) p1_jac (x, t, ret);

  ret.p2.ode = @ p2_ode;

  ret.p2.jac = @ p2_jac;

  ret.p2.sol = @ p2_sol;

  ret.p2.x0 = @ p2_x0;

  function ret = p1_ode (x, t, hook)

    ## Test problem taken from the paper "Multirate Numerical
    ## Integration for Stiff ODEs" by V. Savcenco and R.M.M. Mattheij,
    ## Technische Universiteit Eindhoven, P.O. Box 513, 5600 MB
    ## Eindhoven, The Netherlands. This is a report of December 2008 of
    ## the "Centre for Analysis, Scientific computing and Applications"
    ## (CASA) of this university, currently available at
    ## http://www.win.tue.nl/analysis/reports/rana08-39.pdf. The
    ## Jacobian is non-zero only in the main diagonal and in the
    ## first lower subdiagonal.

    g = hook.p1.helpers.g;

    u = hook.p1.helpers.u;

    ret = 5 - x - 100 * g (cat (1, u (t), x(1 : end - 1)), x);

  endfunction

  function ret = p1_jac (x, t, hook)

    u = hook.p1.helpers.u;

    n = rows (x);

    u_x = x(1 : end - 1);

    ret = - eye (n) - \
	(diag (200 * (max (u_x - 1, 0) - max (u_x - x(2:end) - 1, 0)), \
	       -1) + \
	 diag (200 * max (cat (1, u (t), u_x) - x - 1, 0)));

  endfunction

  function ret = p1_g (u, v)

    ret = max (u - 1, 0) .^ 2 - max (u - v - 1, 0) .^ 2;

  endfunction

  function ret = p1_u (t)

    if (t < 5 || t > 17)
      ret = 0;
    elseif (t <= 10)
      ret = t - 5;
    elseif (t <= 15)
      ret = 5;
    else
      ret = 5 * (17 - t) / 2;
    endif

  endfunction

  function x0 = p1_x0 (n)

    x0 = zeros (2, ceil (n / 2));
    x0(1, :) = 5;
    x0(2, :) = 6.247e-3;
    x0 = x0(:);
    if (rem (n, 2))
      x0 = x0(1 : end - 1);
    endif

  endfunction

  function ret = p2_ode (x, t)

    ## independent differential equations

    ret = - x;

  endfunction

  function ret = p2_jac (x, t)

    ret = - eye (rows (x));

  endfunction

  function ret = p2_sol (x0, t)

    x0 = x0(:);

    t = reshape (t, 1, []);

    a = x0 / exp (- t(1));

    ret = a * exp (- t);

  endfunction

  function x0 = p2_x0 (n)

    x0 = ones (n, 1);

  endfunction

endfunction
