octave_packages/odepkg-0.8.2/odepkg_testsuite_implakzo.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_implakzo (@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 chemical AKZO Nobel testsuite of implicit differential algebraic equations after solving (IDE--test).
  21 %#
  22 %# Run examples with the command
  23 %# @example
  24 %# demo odepkg_testsuite_implakzo
  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_implakzo (vhandle, vrtol)
  33
  34   if (nargin ~= 2) %# Check number and types of all input arguments
  35     help  ('odepkg_testsuite_implakzo');
  36     error ('OdePkg:InvalidArgument', ...
  37            'Number of input arguments must be exactly two');
  38   elseif (~isa (vhandle, 'function_handle') || ~isscalar (vrtol))
  39     print_usage;
  40   end
  41
  42   vret{1} = vhandle; %# The handle for the solver that is used
  43   vret{2} = vrtol;   %# The value for the realtive tolerance
  44   vret{3} = vret{2}; %# The value for the absolute tolerance
  45   vret{4} = vret{2}; %# The value for the first time step
  46   %# Write a debug message on the screen, because this testsuite function
  47   %# may be called more than once from a loop over all solvers present
  48   fprintf (1, ['Testsuite AKZO, testing solver %7s with relative', ...
  49     ' tolerance %2.0e\n'], func2str (vret{1}), vrtol); fflush (1);
  50
  51   %# Setting the integration algorithms option values
  52   vstart = 0.0;   %# The point of time when solving is started
  53   vstop  = 180.0; %# The point of time when solving is stoped
  54   [vinity, vinityd] = odepkg_testsuite_implakzoinit; %# The initial values
  55
  56   vopt = odeset ('Refine', 0, 'RelTol', vret{2}, 'AbsTol', vret{3}, ...
  57     'InitialStep', vret{4}, 'Stats', 'on', 'NormControl', 'off', ...
  58     'Jacobian', @odepkg_testsuite_implakzojac, 'MaxStep', vstop-vstart);
  59     %# ,'OutputFcn', @odeplot, 'MaxStep', 1);
  60
  61   %# Calculate the algorithm, start timer and do solving
  62   tic; vsol = feval (vhandle, @odepkg_testsuite_implakzofun, ...
  63     [vstart, vstop], vinity, vinityd', vopt);
  64   vret{12} = toc;                      %# The value for the elapsed time
  65   vref = odepkg_testsuite_implakzoref; %# Get the reference solution vector
  66   if (exist ('OCTAVE_VERSION') ~= 0)
  67     vlst = vsol.y(end,:);
  68   else
  69     vlst = vsol.y(:,end);
  70   end
  71   vret{5}  = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2});
  72   vret{6}  = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2});
  73   vret{7}  = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals
  74   vret{8}  = vsol.stats.nsteps;   %# The value for success evals
  75   vret{9}  = vsol.stats.nfevals;  %# The value for fun calls
  76   vret{10} = vsol.stats.npds;     %# The value for partial derivations
  77   vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions
  78
  79 %# Return the results for the for the chemical AKZO problem
  80 function res = odepkg_testsuite_implakzofun (t, y, yd, varargin)
  81   k1   = 18.7; k2  = 0.58; k3 = 0.09;   k4  = 0.42;
  82   kbig = 34.4; kla = 3.3;  ks = 115.83; po2 = 0.9;
  83   hen  = 737;
  84
  85   r1  = k1 * y(1)^4 * sqrt (y(2));
  86   r2  = k2 * y(3) * y(4);
  87   r3  = k2 / kbig * y(1) * y(5);
  88   r4  = k3 * y(1) * y(4)^2;
  89   r5  = k4 * y(6)^2 * sqrt (y(2));
  90   fin = kla * (po2 / hen - y(2));
  91
  92   res(1,1) = -2 * r1 + r2 - r3 - r4 - yd(1);
  93   res(2,1) = -0.5 * r1 - r4 - 0.5 * r5 + fin - yd(2);
  94   res(3,1) = r1 - r2 + r3 - yd(3);
  95   res(4,1) = - r2 + r3 - 2 * r4 - yd(4);
  96   res(5,1) = r2 - r3 + r5 - yd(5);
  97   res(6,1) = ks * y(1) * y(4) - y(6) - yd(6);
  98
  99 %# Return the INITIAL values for the chemical AKZO problem
 100 function [y0, yd0] = odepkg_testsuite_implakzoinit ()
 101   y0 = [0.444, 0.00123, 0, 0.007, 0, 115.83 * 0.444 * 0.007];
 102   yd0 = [-0.051, -0.014, 0.025, 0, 0.002, 0];
 103
 104 %# Return the JACOBIAN matrix for the chemical AKZO problem
 105 function [dfdy, dfdyd] = odepkg_testsuite_implakzojac (t, y, varargin)
 106   k1   = 18.7; k2  = 0.58; k3 = 0.09;   k4  = 0.42;
 107   kbig = 34.4; kla = 3.3;  ks = 115.83; po2 = 0.9;
 108   hen  = 737;
 109
 110 %# if (y(2) <= 0)
 111 %#   error ('odepkg_testsuite_implakzojac: Second input argument is negative');
 112 %# end
 113
 114   dfdy = zeros (6, 6);
 115
 116   r11  = 4 * k1 * y(1)^3 * sqrt (y(2));
 117   r12  = 0.5 * k1 * y(1)^4 / sqrt (y(2));
 118   r23  = k2 * y(4);
 119   r24  = k2 * y(3);
 120   r31  = (k2 / kbig) * y(5);
 121   r35  = (k2 / kbig) * y(1);
 122   r41  = k3 * y(4)^2;
 123   r44  = 2 * k3 * y(1) * y(4);
 124   r52  = 0.5 * k4 * y(6)^2 / sqrt (y(2));
 125   r56  = 2 * k4 * y(6) * sqrt (y(2));
 126   fin2 = -kla;
 127
 128   dfdy(1,1) = -2 * r11 - r31 - r41;
 129   dfdy(1,2) = -2 * r12;
 130   dfdy(1,3) = r23;
 131   dfdy(1,4) = r24 - r44;
 132   dfdy(1,5) = -r35;
 133   dfdy(2,1) = -0.5 * r11 - r41;
 134   dfdy(2,2) = -0.5 * r12 - 0.5 * r52 + fin2;
 135   dfdy(2,4) = -r44;
 136   dfdy(2,6) = -0.5 * r56;
 137   dfdy(3,1) = r11 + r31;
 138   dfdy(3,2) = r12;
 139   dfdy(3,3) = -r23;
 140   dfdy(3,4) = -r24;
 141   dfdy(3,5) = r35;
 142   dfdy(4,1) = r31 - 2 * r41;
 143   dfdy(4,3) = -r23;
 144   dfdy(4,4) = -r24 - 2 * r44;
 145   dfdy(4,5) = r35;
 146   dfdy(5,1) = -r31;
 147   dfdy(5,2) = r52;
 148   dfdy(5,3) = r23;
 149   dfdy(5,4) = r24;
 150   dfdy(5,5) = -r35;
 151   dfdy(5,6) = r56;
 152   dfdy(6,1) = ks * y(4);
 153   dfdy(6,4) = ks * y(1);
 154   dfdy(6,6) = -1;
 155
 156   dfdyd = - [ 1, 0, 0, 0, 0, 0;
 157               0, 1, 0, 0, 0, 0;
 158               0, 0, 1, 0, 0, 0;
 159               0, 0, 0, 1, 0, 0;
 160               0, 0, 0, 0, 1, 0;
 161               0, 0, 0, 0, 0, 1 ];
 162
 163 %# For the implicit form of the chemical AKZO Nobel problem a mass
 164 %# matrix is not needed. This mass matrix is needed if the problem
 165 %# is formulated in explicit form (cf. odepkg_testsuite_cemakzo.m).
 166 %# function mass = odepkg_testsuite_implakzomass (t, y, varargin)
 167 %#   mass =  [ 1, 0, 0, 0, 0, 0;
 168 %#             0, 1, 0, 0, 0, 0;
 169 %#             0, 0, 1, 0, 0, 0;
 170 %#             0, 0, 0, 1, 0, 0;
 171 %#             0, 0, 0, 0, 1, 0;
 172 %#             0, 0, 0, 0, 0, 0 ];
 173
 174 %# Return the REFERENCE values for the chemical AKZO problem
 175 function y = odepkg_testsuite_implakzoref ()
 176   y(1,1) = 0.11507949206617e+0;
 177   y(2,1) = 0.12038314715677e-2;
 178   y(3,1) = 0.16115628874079e+0;
 179   y(4,1) = 0.36561564212492e-3;
 180   y(5,1) = 0.17080108852644e-1;
 181   y(6,1) = 0.48735313103074e-2;
 182
 183 %!demo
 184 %! vsolver = {@odebdi};
 185 %! for vcnt=1:length (vsolver)
 186 %!   vakzo{vcnt,1} = odepkg_testsuite_implakzo (vsolver{vcnt}, 1e-7);
 187 %! end
 188 %! vakzo
 189
 190 %# Local Variables: ***
 191 %# mode: octave ***
 192 %# End: ***