octave_packages/odepkg-0.8.2/odepkg_testsuite_implrober.m

   1 %# Copyright (C) 2007-2012, Thomas Treichl <treichl@users.sourceforge.net>
   2 %# OdePkg - A package for solving ordinary differential equations and more
   3 %#
   4 %# This program is free software; you can redistribute it and/or modify
   5 %# it under the terms of the GNU General Public License as published by
   6 %# the Free Software Foundation; either version 2 of the License, or
   7 %# (at your option) any later version.
   8 %#
   9 %# This program is distributed in the hope that it will be useful,
  10 %# but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 %# GNU General Public License for more details.
  13 %#
  14 %# You should have received a copy of the GNU General Public License
  15 %# along with this program; If not, see <http://www.gnu.org/licenses/>.
  16
  17 %# -*- texinfo -*-
  18 %# @deftypefn {Function File} {[@var{solution}] =} odepkg_testsuite_implrober (@var{@@solver}, @var{reltol})
  19 %#
  20 %# If this function is called with two input arguments and the first input argument @var{@@solver} is a function handle describing an OdePkg solver and the second input argument @var{reltol} is a double scalar describing the relative error tolerance then return a cell array @var{solution} with performance informations about the implicit form of the modified ROBERTSON testsuite of implicit differential algebraic equations after solving (IDE--test).
  21 %#
  22 %# Run examples with the command
  23 %# @example
  24 %# demo odepkg_testsuite_implrober
  25 %# @end example
  26 %#
  27 %# This function has been ported from the "Test Set for IVP solvers" which is developed by the INdAM Bari unit project group "Codes and Test Problems for Differential Equations", coordinator F. Mazzia.
  28 %# @end deftypefn
  29 %#
  30 %# @seealso{odepkg}
  31
  32 function vret = odepkg_testsuite_implrober (vhandle, vrtol)
  33
  34   %# Check number and types of all input arguments
  35   if (nargin ~= 2)
  36     help  ('odepkg_testsuite_implrober');
  37     error ('OdePkg:InvalidArgument', ...
  38            'Number of input arguments must be exactly two');
  39   elseif (~isa (vhandle, 'function_handle') || ~isscalar (vrtol))
  40     print_usage;
  41   end
  42
  43   vret{1} = vhandle; %# The handle for the solver that is used
  44   vret{2} = vrtol;   %# The value for the realtive tolerance
  45   vret{3} = vret{2} * 1e-2; %# The value for the absolute tolerance
  46   vret{4} = vret{2}; %# The value for the first time step
  47   %# Write a debug message on the screen, because this testsuite function
  48   %# may be called more than once from a loop over all present solvers
  49   fprintf (1, ['Testsuite implicit ROBERTSON, testing solver %7s with relative', ...
  50      ' tolerance %2.0e\n'], func2str (vret{1}), vrtol); fflush (1);
  51
  52   %# Setting the integration algorithm options
  53   vstart = 0;    %# The point of time when solving is started
  54   vstop  = 1e11; %# The point of time when solving is stoped
  55   [vinity, vinityd] = odepkg_testsuite_implroberinit; %# The initial values
  56
  57   vopt = odeset ('Refine', 0, 'RelTol', vret{2}, 'AbsTol', vret{3}, ...
  58     'InitialStep', vret{4}, 'Stats', 'on', 'NormControl', 'off', ...
  59     'Jacobian', @odepkg_testsuite_implroberjac, 'MaxStep', vstop-vstart);
  60     %# 'OutputFcn', @odeplot);
  61
  62   %# Calculate the algorithm, start timer and do solving
  63   tic; vsol = feval (vhandle, @odepkg_testsuite_implroberfun, ...
  64     [vstart, vstop], vinity, vinityd', vopt);
  65   vret{12} = toc;                       %# The value for the elapsed time
  66   vref = odepkg_testsuite_implroberref; %# Get the reference solution vector
  67   if (exist ('OCTAVE_VERSION') ~= 0)
  68     vlst = vsol.y(end,:);
  69   else
  70     vlst = vsol.y(:,end);
  71   end
  72   vret{5}  = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2});
  73   vret{6}  = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2});
  74   vret{7}  = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals
  75   vret{8}  = vsol.stats.nsteps;   %# The value for success evals
  76   vret{9}  = vsol.stats.nfevals;  %# The value for fun calls
  77   vret{10} = vsol.stats.npds;     %# The value for partial derivations
  78   vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions
  79
  80 %#function odepkg_testsuite_implrober ()
  81 %#  A = odeset ('RelTol', 1e-4, ... %# proprietary ode15i needs 1e-6 to be stable
  82 %#             'AbsTol', [1e-6, 1e-10, 1e-6], ...
  83 %#             'Jacobian', @odepkg_testsuite_implroberjac);
  84 %#  [y0, yd0] = odepkg_testsuite_implroberinit;
  85 %#  odebdi (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0', A)
  86 %#  [y0, yd0] = odepkg_testsuite_implroberinit;
  87 %#  odebdi (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0')
  88
  89 %# Return the results for the for the implicit ROBERTSON problem
  90 function res = odepkg_testsuite_implroberfun (t, y, yd, varargin)
  91   res(1,1) = -0.04 * y(1) + 1e4 * y(2) * y(3) - yd(1);
  92   res(2,1) =  0.04 * y(1) - 1e4 * y(2) * y(3) - 3e7 * y(2)^2 - yd(2);
  93   res(3,1) =  y(1) + y(2) + y(3) - 1;
  94
  95 %# Return the INITIAL values for the implicit ROBERTSON problem
  96 function [y0, yd0] = odepkg_testsuite_implroberinit ()
  97   y0 = [1, 0, 0];
  98   yd0 = [-4e-2, 4e-2, 0];
  99
 100 %# Return the JACOBIAN matrix for the implicit ROBERTSON problem
 101 function [dfdy, dfdyd] = odepkg_testsuite_implroberjac (t, y, yd, varargin)
 102   dfdy(1,1)  = -0.04;
 103   dfdy(1,2)  =  1e4 * y(3);
 104   dfdy(1,3)  =  1e4 * y(2);
 105   dfdy(2,1)  =  0.04;
 106   dfdy(2,2)  = -1e4 * y(3) - 6e7 * y(2);
 107   dfdy(2,3)  = -1e4 * y(2);
 108   dfdy(3,1)  =  1;
 109   dfdy(3,2)  =  1;
 110   dfdy(3,3)  =  1;
 111
 112   dfdyd(1,1) = -1;
 113   dfdyd(2,2) = -1;
 114   dfdyd(3,3) =  0;
 115
 116 %# For the implicit form of the Robertson problem a mass matrix is not
 117 %# allowed. This mass matrix is only needed if the Robertson problem
 118 %# is formulated in explicit form (cf. odepkg_testsuite_implrober.m).
 119 %# function mass = odepkg_testsuite_implrobermass (t, y, varargin)
 120 %#   mass =  [1, 0, 0; 0, 1, 0; 0, 0, 0];
 121
 122 %# Return the REFERENCE values for the implicit ROBERTSON problem
 123 function y = odepkg_testsuite_implroberref ()
 124   y(1) =  0.20833401497012e-07;
 125   y(2) =  0.83333607703347e-13;
 126   y(3) =  0.99999997916650e+00;
 127
 128 %!demo
 129 %! vsolver = {@odebdi};
 130 %! for vcnt=1:length (vsolver)
 131 %!   virob{vcnt,1} = odepkg_testsuite_implrober (vsolver{vcnt}, 1e-7);
 132 %! end
 133 %! virob
 134
 135 %# Local Variables: ***
 136 %# mode: octave ***
 137 %# End: ***