octave_packages/odepkg-0.8.2/odepkg_testsuite_pollution.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_pollution (@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 the cell array @var{solution} with performance informations about the POLLUTION testsuite of ordinary differential equations after solving (ODE--test).
  21 %#
  22 %# Run examples with the command
  23 %# @example
  24 %# demo odepkg_testsuite_pollution
  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_pollution (vhandle, vrtol)
  33
  34   if (nargin ~= 2) %# Check number and types of all input arguments
  35     help  ('odepkg_testsuite_pollution');
  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 POLLUTION, 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  = 60.0;  %# The point of time when solving is stoped
  54   vinit  = odepkg_testsuite_pollutioninit; %# 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_pollutionjac, 'MaxStep', vstop-vstart);
  59
  60   %# Calculate the algorithm, start timer and do solving
  61   tic; vsol = feval (vhandle, @odepkg_testsuite_pollutionfun, ...
  62     [vstart, vstop], vinit, vopt);
  63   vret{12} = toc;                       %# The value for the elapsed time
  64   vref = odepkg_testsuite_pollutionref; %# Get the reference solution vector
  65   if (exist ('OCTAVE_VERSION') ~= 0)
  66     vlst = vsol.y(end,:);
  67   else
  68     vlst = vsol.y(:,end);
  69   end
  70   vret{5}  = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2});
  71   vret{6}  = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2});
  72   vret{7}  = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals
  73   vret{8}  = vsol.stats.nsteps;   %# The value for success evals
  74   vret{9}  = vsol.stats.nfevals;  %# The value for fun calls
  75   vret{10} = vsol.stats.npds;     %# The value for partial derivations
  76   vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions
  77
  78 %# Returns the results for the for the POLLUTION problem
  79 function f = odepkg_testsuite_pollutionfun (t, y, varargin)
  80   f(01,1) = - 0.350 * y(1) + 0.266e2 * y(2) * y(4) + 0.123e5 * y(2) * y(5) + ...
  81               0.165e5 * y(2) * y(11) - 0.900e4 * y(1) * y(11) + 0.220e-1 * y(13) + ...
  82               0.120e5 * y(2) * y(10) - 0.163e5 * y(1) * y(6) + 0.578e1 * y(19) - ...
  83               0.474e-1 * y(1) * y(4) - 0.178e4 * y(1) * y(19) + 0.312e1 * y(20);
  84   f(02,1) = + 0.350 * y(1) - 0.266e2 * y(2) * y(4) - 0.123e5 * y(2) * y(5) - ...
  85               0.165e5 * y(2) * y(11) - 0.120e5 * y(2) * y(10) + 0.210e1 * y(19);
  86   f(03,1) = + 0.350 * y(1) - 0.480e7 * y(3) + 0.175e-1 * y(4) + 0.444e12 * y(16) + 0.578e1 * y(19);
  87   f(04,1) = - 0.266e2 * y(2) * y(4) + 0.480e7 * y(3) - 0.350e-3 * y(4) - ...
  88               0.175e-1 * y(4) - 0.474e-1 * y(1) * y(4);
  89   f(05,1) = - 0.123e5 * y(2) * y(5) + 2*0.860e-3 * y(7) + 0.150e5 * y(6) * y(7) + ...
  90               0.130e-3 * y(9) + 0.188e1 * y(14) + 0.124e4 * y(6) * y(17);
  91   f(06,1) = + 0.123e5 * y(2) * y(5) - 0.150e5 * y(6) * y(7) - 0.240e5 * y(6) * y(9) - ...
  92               0.163e5 * y(1) * y(6) + 2*0.100e9 * y(16) - 0.124e4 * y(6) * y(17);
  93   f(07,1) = - 0.860e-3 * y(7) - 0.820e-3 * y(7) - 0.150e5 * y(6) * y(7) + 0.188e1 * y(14);
  94   f(08,1) = + 0.860e-3 * y(7) + 0.820e-3 * y(7) + 0.150e5 * y(6) * y(7) + 0.130e-3 * y(9);
  95   f(09,1) = - 0.130e-3 * y(9) - 0.240e5 * y(6) * y(9);
  96   f(10,1) = + 0.130e-3 * y(9) + 0.165e5 * y(2) * y(11) - 0.120e5 * y(2) * y(10);
  97   f(11,1) = + 0.240e5 * y(6) * y(9) - 0.165e5 * y(2) * y(11) - 0.900e4 * y(1) * y(11) + ...
  98               0.220e-1 * y(13);
  99   f(12,1) = + 0.165e5 * y(2) * y(11);
 100   f(13,1) = + 0.900e4 * y(1) * y(11) - 0.220e-1 * y(13);
 101   f(14,1) = + 0.120e5 * y(2) * y(10) - 0.188e1 * y(14);
 102   f(15,1) = + 0.163e5 * y(1) * y(6);
 103   f(16,1) = + 0.350e-3 * y(4) - 0.100e9 * y(16) - 0.444e12 * y(16);
 104   f(17,1) = - 0.124e4 * y(6) * y(17);
 105   f(18,1) = + 0.124e4 * y(6) * y(17);
 106   f(19,1) = - 0.210e1 * y(19) - 0.578e1 * y(19) + 0.474e-1 * y(1) * y(4) - ...
 107               0.178e4 * y(1) * y(19) + 0.312e1 * y(20);
 108   f(20,1) = + 0.178e4 * y(1) * y(19) - 0.312e1 * y(20);
 109
 110 %# Returns the INITIAL values for the POLLUTION problem
 111 function vinit = odepkg_testsuite_pollutioninit ()
 112   vinit = [0, 0.2, 0, 0.04, 0, 0, 0.1, 0.3, 0.01, ...
 113            0, 0, 0, 0, 0, 0, 0, 0.007, 0, 0, 0];
 114
 115 %# Returns the JACOBIAN matrix for the POLLUTION problem
 116 function dfdy = odepkg_testsuite_pollutionjac (t, y)
 117   k1  = 0.35e0;   k2  = 0.266e2; k3  = 0.123e5;  k4  = 0.86e-3;
 118   k5  = 0.82e-3;  k6  = 0.15e5;  k7  = 0.13e-3;  k8  = 0.24e5;
 119   k9  = 0.165e5;  k10 = 0.9e4;   k11 = 0.22e-1;  k12 = 0.12e5;
 120   k13 = 0.188e1;  k14 = 0.163e5; k15 = 0.48e7;   k16 = 0.35e-3;
 121   k17 = 0.175e-1; k18 = 0.1e9;   k19 = 0.444e12; k20 = 0.124e4;
 122   k21 = 0.21e1;   k22 = 0.578e1; k23 = 0.474e-1; k24 = 0.178e4;
 123   k25 = 0.312e1;
 124
 125   dfdy(1,1)   = -k1 - k10 * y(11) - k14 * y(6) - k23 * y(4) - k24 * y(19);
 126   dfdy(1,11)  = -k10 * y(1) + k9 * y(2);
 127   dfdy(1,6)   = -k14 * y(1);
 128   dfdy(1,4)   = -k23 * y(1) + k2 * y(2);
 129   dfdy(1,19)  = -k24 * y(1) + k22;
 130   dfdy(1,2)   =  k2 * y(4) + k9 * y(11) + k3 * y(5) + k12 * y(10);
 131   dfdy(1,13)  =  k11;
 132   dfdy(1,20)  =  k25;
 133   dfdy(1,5)   =  k3 * y(2);
 134   dfdy(1,10)  =  k12 * y(2);
 135
 136   dfdy(2,4)   = -k2 * y(2);
 137   dfdy(2,5)   = -k3 * y(2);
 138   dfdy(2,11)  = -k9 * y(2);
 139   dfdy(2,10)  = -k12 * y(2);
 140   dfdy(2,19)  =  k21;
 141   dfdy(2,1)   =  k1;
 142   dfdy(2,2)   = -k2 * y(4) - k3 * y(5) - k9 * y(11) - k12 * y(10);
 143
 144   dfdy(3,1)   =  k1;
 145   dfdy(3,4)   =  k17;
 146   dfdy(3,16)  =  k19;
 147   dfdy(3,19)  =  k22;
 148   dfdy(3,3)   = -k15;
 149
 150   dfdy(4,4)   = -k2 * y(2) - k16 - k17 - k23 * y(1);
 151   dfdy(4,2)   = -k2 * y(4);
 152   dfdy(4,1)   = -k23 * y(4);
 153   dfdy(4,3)   =  k15;
 154
 155   dfdy(5,5)   = -k3 * y(2);
 156   dfdy(5,2)   = -k3 * y(5);
 157   dfdy(5,7)   =  2d0 * k4 + k6 * y(6);
 158   dfdy(5,6)   =  k6 * y(7) + k20 * y(17);
 159   dfdy(5,9)   =  k7;
 160   dfdy(5,14)  =  k13;
 161   dfdy(5,17)  =  k20 * y(6);
 162
 163   dfdy(6,6)   = -k6 * y(7) - k8 * y(9) - k14 * y(1) - k20 * y(17);
 164   dfdy(6,7)   = -k6 * y(6);
 165   dfdy(6,9)   = -k8 * y(6);
 166   dfdy(6,1)   = -k14 * y(6);
 167   dfdy(6,17)  = -k20 * y(6);
 168   dfdy(6,2)   =  k3 * y(5);
 169   dfdy(6,5)   =  k3 * y(2);
 170   dfdy(6,16)  =  2d0 * k18;
 171
 172   dfdy(7,7)   = -k4 - k5 - k6 * y(6);
 173   dfdy(7,6)   = -k6 * y(7);
 174   dfdy(7,14)  =  k13;
 175
 176   dfdy(8,7)   =  k4 + k5 + k6 * y(6);
 177   dfdy(8,6)   =  k6 * y(7);
 178   dfdy(8,9)   =  k7;
 179
 180   dfdy(9,9)   = -k7 - k8 * y(6);
 181   dfdy(9,6)   = -k8 * y(9);
 182
 183   dfdy(10,10) = -k12 * y(2);
 184   dfdy(10,2)  = -k12 * y(10) + k9 * y(11);
 185   dfdy(10,9)  =  k7;
 186   dfdy(10,11) =  k9 * y(2);
 187
 188   dfdy(11,11) = -k9 * y(2) - k10 * y(1);
 189   dfdy(11,2)  = -k9 * y(11);
 190   dfdy(11,1)  = -k10 * y(11);
 191   dfdy(11,9)  =  k8 * y(6);
 192   dfdy(11,6)  =  k8 * y(9);
 193   dfdy(11,13) =  k11;
 194
 195   dfdy(12,11) =  k9 * y(2);
 196   dfdy(12,2)  =  k9 * y(11);
 197
 198   dfdy(13,13) = -k11;
 199   dfdy(13,11) =  k10 * y(1);
 200   dfdy(13,1)  =  k10 * y(11);
 201
 202   dfdy(14,14) = -k13;
 203   dfdy(14,10) =  k12 * y(2);
 204   dfdy(14,2)  =  k12 * y(10);
 205
 206   dfdy(15,1)  =  k14 * y(6);
 207   dfdy(15,6)  =  k14 * y(1);
 208
 209   dfdy(16,16) = -k18 - k19;
 210   dfdy(16,4)  =  k16;
 211
 212   dfdy(17,17) = -k20 * y(6);
 213   dfdy(17,6)  = -k20 * y(17);
 214
 215   dfdy(18,17) =  k20 * y(6);
 216   dfdy(18,6)  =  k20 * y(17);
 217
 218   dfdy(19,19) = -k21 - k22 - k24 * y(1);
 219   dfdy(19,1)  = -k24 * y(19) + k23 * y(4);
 220   dfdy(19,4)  =  k23 * y(1);
 221   dfdy(19,20) =  k25;
 222
 223   dfdy(20,20) = -k25;
 224   dfdy(20,1)  =  k24 * y(19);
 225   dfdy(20,19) =  k24 * y(1);
 226
 227 %# Returns the REFERENCE values for the POLLUTION problem
 228 function y = odepkg_testsuite_pollutionref ()
 229   y(01,1) = 0.56462554800227 * 10^(-1);
 230   y(02,1) = 0.13424841304223 * 10^(+0);
 231   y(03,1) = 0.41397343310994 * 10^(-8);
 232   y(04,1) = 0.55231402074843 * 10^(-2);
 233   y(05,1) = 0.20189772623021 * 10^(-6);
 234   y(06,1) = 0.20189772623021 * 10^(-6);
 235   y(07,1) = 0.77842491189979 * 10^(-1);
 236   y(08,1) = 0.77842491189979 * 10^(+0);
 237   y(09,1) = 0.74940133838804 * 10^(-2);
 238   y(10,1) = 0.16222931573015 * 10^(-7);
 239   y(11,1) = 0.11358638332570 * 10^(-7);
 240   y(12,1) = 0.22305059757213 * 10^(-2);
 241   y(13,1) = 0.20871628827986 * 10^(-3);
 242   y(14,1) = 0.13969210168401 * 10^(-4);
 243   y(15,1) = 0.89648848568982 * 10^(-2);
 244   y(16,1) = 0.43528463693301 * 10^(-17);
 245   y(17,1) = 0.68992196962634 * 10^(-2);
 246   y(18,1) = 0.10078030373659 * 10^(-3);
 247   y(19,1) = 0.17721465139699 * 10^(-5);
 248   y(20,1) = 0.56829432923163 * 10^(-4);
 249
 250 %!demo
 251 %! %% vsolver = {@ode23, @ode45, @ode54, @ode78, ...
 252 %! %%   @odebda, @oders, @ode2r, @ode5r, @odesx};
 253 %! vsolver = {@odebda, @oders, @ode2r, @ode5r, @odesx};
 254 %! for vcnt=1:length (vsolver)
 255 %!   poll{vcnt,1} = odepkg_testsuite_pollution (vsolver{vcnt}, 1e-7);
 256 %! end
 257 %! poll
 258
 259 %# Local Variables: ***
 260 %# mode: octave ***
 261 %# End: ***