Mercurial > forge
changeset 4012:d573d276bc63 octave-forge
More testsuite functions for IDE solvers.
| author | treichl |
|---|---|
| date | Fri, 09 Nov 2007 16:35:11 +0000 |
| parents | 19755d677c22 |
| children | ab34413d4486 |
| files | main/odepkg/inst/odepkg_testsuite_implakzo.m main/odepkg/inst/odepkg_testsuite_implrober.m |
| diffstat | 2 files changed, 333 insertions(+), 143 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/odepkg/inst/odepkg_testsuite_implakzo.m Fri Nov 09 16:35:11 2007 +0000 @@ -0,0 +1,194 @@ +%# Copyright (C) 2007, Thomas Treichl <treichl@users.sourceforge.net> +%# OdePkg - Package for solving ordinary differential equations with octave +%# +%# 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 2 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, write to the Free Software +%# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + +%# -*- texinfo -*- +%# @deftypefn {Function File} {[@var{solution}] =} odepkg_testsuite_implakzo (@var{@@solver}, @var{reltol}) +%# +%# 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). +%# +%# Run examples with the command +%# @example +%# demo odepkg_testsuite_implakzo +%# @end example +%# @end deftypefn +%# +%# @seealso{odepkg} + +function vret = odepkg_testsuite_implakzo (vhandle, vrtol) + + if (nargin ~= 2) %# Check number and types of all input arguments + help ('odepkg_testsuite_implakzo'); + error ('OdePkg:InvalidInputArgument', ... + 'Number of input arguments must be exactly two'); + elseif (~isa (vhandle, 'function_handle') || ~isscalar (vrtol)) + print_usage; + end + + vret{1} = vhandle; %# The handle for the solver that is used + vret{2} = vrtol; %# The value for the realtive tolerance + vret{3} = vret{2}; %# The value for the absolute tolerance + vret{4} = vret{2}; %# The value for the first time step + %# Write a debug message on the screen, because this testsuite function + %# may be called more than once from a loop over all solvers present + fprintf (1, ['Testsuite AKZO, testing solver %7s with relative', ... + ' tolerance %2.0e\n'], func2str (vret{1}), vrtol); fflush (1); + + %# Setting the integration algorithms option values + vstart = 0.0; %# The point of time when solving is started + vstop = 180.0; %# The point of time when solving is stoped + [vinity, vinityd] = odepkg_testsuite_implakzoinit; %# The initial values + + vopt = odeset ('Refine', 0, 'RelTol', vret{2}, 'AbsTol', vret{3}, ... + 'InitialStep', vret{4}, 'Stats', 'on', 'NormControl', 'off', ... + 'Jacobian', @odepkg_testsuite_implakzojac, 'MaxStep', vstop-vstart); + %# 'OutputFcn', @odeplot); + + %# Calculate the algorithm, start timer and do solving + tic; vsol = feval (vhandle, @odepkg_testsuite_implakzofun, ... + [vstart, vstop], vinity, vinityd', vopt); + vret{12} = toc; %# The value for the elapsed time + vref = odepkg_testsuite_implakzoref; %# Get the reference solution vector + if (max (size (vsol.y(end,:))) == max (size (vref))), vlst = vsol.y(end,:); + elseif (max (size (vsol.y(:,end))) == max (size (vref))), vlst = vsol.y(:,end); + end + vret{5} = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2}); + vret{6} = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2}); +%# vret{7} = vsol.stats.success + vsol.stats.failed; %# The value for all evals +%# vret{8} = vsol.stats.success; %# The value for success evals +%# vret{9} = vsol.stats.fevals; %# The value for fun calls +%# vret{10} = vsol.stats.partial; %# The value for partial derivations +%# vret{11} = vsol.stats.ludecom; %# The value for LU decompositions + vret{7} = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals + vret{8} = vsol.stats.nsteps; %# The value for success evals + vret{9} = vsol.stats.nfevals; %# The value for fun calls + vret{10} = vsol.stats.nsolves; %# The value for partial derivations + vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions + +%# Return the results for the for the chemical AKZO problem +function res = odepkg_testsuite_implakzofun (t, y, yd, varargin) + k1 = 18.7; k2 = 0.58; k3 = 0.09; k4 = 0.42; + kbig = 34.4; kla = 3.3; ks = 115.83; po2 = 0.9; + hen = 737; + + r1 = k1 * y(1)^4 * sqrt (y(2)); + r2 = k2 * y(3) * y(4); + r3 = k2 / kbig * y(1) * y(5); + r4 = k3 * y(1) * y(4)^2; + r5 = k4 * y(6)^2 * sqrt (y(2)); + fin = kla * (po2 / hen - y(2)); + + res(1,1) = -2 * r1 + r2 - r3 - r4 - yd(1); + res(2,1) = -0.5 * r1 - r4 - 0.5 * r5 + fin - yd(2); + res(3,1) = r1 - r2 + r3 - yd(3); + res(4,1) = - r2 + r3 - 2 * r4 - yd(4); + res(5,1) = r2 - r3 + r5 - yd(5); + res(6,1) = ks * y(1) * y(4) - y(6) - yd(6); + +%# Return the INITIAL values for the chemical AKZO problem +function [y0, yd0] = odepkg_testsuite_implakzoinit () + y0 = [0.444, 0.00123, 0, 0.007, 0, 115.83 * 0.444 * 0.007]; + yd0 = [-0.051, -0.014, 0.025, 0, 0.002, 0]; + +%# Return the JACOBIAN matrix for the chemical AKZO problem +function [dfdy, dfdyd] = odepkg_testsuite_implakzojac (t, y, varargin) + k1 = 18.7; k2 = 0.58; k3 = 0.09; k4 = 0.42; + kbig = 34.4; kla = 3.3; ks = 115.83; po2 = 0.9; + hen = 737; + +%# if (y(2) <= 0) +%# error ('odepkg_testsuite_implakzojac: Second input argument is negative'); +%# end + + dfdy = zeros (6, 6); + + r11 = 4 * k1 * y(1)^3 * sqrt (y(2)); + r12 = 0.5 * k1 * y(1)^4 / sqrt (y(2)); + r23 = k2 * y(4); + r24 = k2 * y(3); + r31 = (k2 / kbig) * y(5); + r35 = (k2 / kbig) * y(1); + r41 = k3 * y(4)^2; + r44 = 2 * k3 * y(1) * y(4); + r52 = 0.5 * k4 * y(6)^2 / sqrt (y(2)); + r56 = 2 * k4 * y(6) * sqrt (y(2)); + fin2 = -kla; + + dfdy(1,1) = -2 * r11 - r31 - r41; + dfdy(1,2) = -2 * r12; + dfdy(1,3) = r23; + dfdy(1,4) = r24 - r44; + dfdy(1,5) = -r35; + dfdy(2,1) = -0.5 * r11 - r41; + dfdy(2,2) = -0.5 * r12 - 0.5 * r52 + fin2; + dfdy(2,4) = -r44; + dfdy(2,6) = -0.5 * r56; + dfdy(3,1) = r11 + r31; + dfdy(3,2) = r12; + dfdy(3,3) = -r23; + dfdy(3,4) = -r24; + dfdy(3,5) = r35; + dfdy(4,1) = r31 - 2 * r41; + dfdy(4,3) = -r23; + dfdy(4,4) = -r24 - 2 * r44; + dfdy(4,5) = r35; + dfdy(5,1) = -r31; + dfdy(5,2) = r52; + dfdy(5,3) = r23; + dfdy(5,4) = r24; + dfdy(5,5) = -r35; + dfdy(5,6) = r56; + dfdy(6,1) = ks * y(4); + dfdy(6,4) = ks * y(1); + dfdy(6,6) = -1; + + dfdyd = - [ 1, 0, 0, 0, 0, 0; + 0, 1, 0, 0, 0, 0; + 0, 0, 1, 0, 0, 0; + 0, 0, 0, 1, 0, 0; + 0, 0, 0, 0, 1, 0; + 0, 0, 0, 0, 0, 1 ]; + +%# For the implicit form of the chemical AKZO Nobel problem a mass +%# matrix is not needed. This mass matrix is needed if the problem +%# is formulated in explicit form (cf. odepkg_testsuite_cemakzo.m). +%# function mass = odepkg_testsuite_implakzomass (t, y, varargin) +%# mass = [ 1, 0, 0, 0, 0, 0; +%# 0, 1, 0, 0, 0, 0; +%# 0, 0, 1, 0, 0, 0; +%# 0, 0, 0, 1, 0, 0; +%# 0, 0, 0, 0, 1, 0; +%# 0, 0, 0, 0, 0, 0 ]; + +%# Return the REFERENCE values for the chemical AKZO problem +function y = odepkg_testsuite_implakzoref () + y(1,1) = 0.11507949206617e+0; + y(2,1) = 0.12038314715677e-2; + y(3,1) = 0.16115628874079e+0; + y(4,1) = 0.36561564212492e-3; + y(5,1) = 0.17080108852644e-1; + y(6,1) = 0.48735313103074e-2; + +%!demo +%! vsolver = {@odebdi}; +%! for vcnt=1:length (vsolver) +%! vakzo{vcnt,1} = odepkg_testsuite_implakzo (vsolver{vcnt}, 1e-7); +%! end +%! vakzo + +%# Local Variables: *** +%# mode: octave *** +%# End: ***
--- a/main/odepkg/inst/odepkg_testsuite_implrober.m Fri Nov 09 16:18:33 2007 +0000 +++ b/main/odepkg/inst/odepkg_testsuite_implrober.m Fri Nov 09 16:35:11 2007 +0000 @@ -1,143 +1,139 @@ -%# Copyright (C) 2006, Thomas Treichl <treichl@users.sourceforge.net> -%# OdePkg - Package for solving ordinary differential equations with octave -%# -%# 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 2 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, write to the Free Software -%# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA - -%# -*- texinfo -*- -%# @deftypefn {Function File} {[@var{solution}] =} odepkg_testsuite_implrober (@var{@@solver}, @var{reltol}) -%# -%# 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). -%# -%# Run examples with the command -%# @example -%# demo odepkg_testsuite_implrober -%# @end example -%# @end deftypefn -%# -%# @seealso{odepkg} - -function vret = odepkg_testsuite_implrober (vhandle, vrtol) - - %# Check number and types of all input arguments - if (nargin ~= 2) - help ('odepkg_testsuite_implrober'); - error ('OdePkg:odepkg_testsuite_implrober:InvalidInputArgument', ... - 'Number of input arguments must be exactly two'); - elseif (~isa (vhandle, 'function_handle') || ~isscalar (vrtol)) - usage ('odepkg_testsuite_implrober (@solver, reltol)'); - end - - vret{1} = vhandle; %# The name for the solver that is used - vret{2} = vrtol; %# The value for the realtive tolerance - vret{3} = vret{2} * 1e-2; %# The value for the absolute tolerance - vret{4} = vret{2}; %# The value for the first time step - %# Write a debug message on the screen, because this testsuite function - %# may be called more than once from a loop over all present solvers - fprintf (1, ['Testsuite implicit ROBERTSON, testing solver %7s with relative', ... - ' tolerance %2.0e\n'], func2str (vret{1}), vrtol); fflush (1); - - %# Setting the integration algorithm options - vstart = 0; %# The point of time when solving is started - vstop = 1e11; %# The point of time when solving is stoped - [vinity, vinityd] = odepkg_testsuite_implroberinit; %# The initial values - %# Cf. note for mass matrix below at commented out function definition - %# vmass = odepkg_testsuite_implrobermass; %# The mass matrix - - vopt = odeset ('Refine', 0, 'RelTol', vret{2}, 'AbsTol', vret{3}, ... - 'InitialStep', vret{4}, 'Stats', 'on', 'NormControl', 'off', ... - 'Jacobian', @odepkg_testsuite_implroberjac, ... %# 'Mass', vmass - 'MaxStep', vstop-vstart); - - %# Calculate the algorithm, start timer and do solving - tic; vsol = feval (vhandle, @odepkg_testsuite_implroberfun, ... - [vstart, vstop], vinity, vinityd', vopt); - vsol.x(end) - vsol.y(:,end) - vret{12} = toc; %# The value for the elapsed time - vref = odepkg_testsuite_implroberref; %# Get the reference solution vector - if (max (size (vsol.y(end,:))) == max (size (vref))), vlst = vsol.y(end,:); - elseif (max (size (vsol.y(:,end))) == max (size (vref))), vlst = vsol.y(:,end); - end - vret{5} = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2}); - vret{6} = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2}); - %# vret{7} = vsol.stats.success + vsol.stats.failed; %# The value for all evals - %# vret{8} = vsol.stats.success; %# The value for success evals - %# vret{9} = vsol.stats.fevals; %# The value for fun calls - %# vret{10} = vsol.stats.partial; %# The value for partial derivations - %# vret{11} = vsol.stats.ludecom; %# The value for LU decompositions - vret{7} = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals - vret{8} = vsol.stats.nsteps; %# The value for success evals - vret{9} = vsol.stats.nfevals; %# The value for fun calls - vret{10} = vsol.stats.nsolves; %# The value for partial derivations - vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions - -%#function odepkg_testsuite_implrober () -%# A = odeset ('RelTol', 1e-4, ... %# MATLAB ode15i needs 1e-6 to be stable -%# 'AbsTol', [1e-6, 1e-10, 1e-6], ... -%# 'Jacobian', @odepkg_testsuite_implroberjac); -%# [y0, yd0] = odepkg_testsuite_implroberinit; -%# ode15i (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0', A) -%# [y0, yd0] = odepkg_testsuite_implroberinit; -%# ode15i (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0') - -%# Returns the results for the for the implicit ROBERTSON problem -function res = odepkg_testsuite_implroberfun (t, y, yd, varargin) - res(1,1) = -0.04 * y(1) + 1e4 * y(2) * y(3) - yd(1); - res(2,1) = 0.04 * y(1) - 1e4 * y(2) * y(3) - 3e7 * y(2)^2 - yd(2); - res(3,1) = y(1) + y(2) + y(3) - 1; - -%# Returns the INITIAL values for the implicit ROBERTSON problem -function [y0, yd0] = odepkg_testsuite_implroberinit () - y0 = [1, 0, 0]; - yd0 = [-4e-2, 4e-2, 0]; - -%# Returns the JACOBIAN matrix for the implicit ROBERTSON problem -function [dfdy, dfdyd] = odepkg_testsuite_implroberjac (t, y, yd, varargin) - dfdy(1,1) = -0.04; - dfdy(1,2) = 1e4 * y(3); - dfdy(1,3) = 1e4 * y(2); - dfdy(2,1) = 0.04; - dfdy(2,2) = -1e4 * y(3) - 6e7 * y(2); - dfdy(2,3) = -1e4 * y(2); - dfdy(3,1) = 1; - dfdy(3,2) = 1; - dfdy(3,3) = 1; - - dfdyd(1,1) = -1; - dfdyd(2,2) = -1; - dfdyd(3,3) = 0; - -%# For the implicit form of the Robertson problem a mass matrix is not -%# allowed. This mass matrix is only needed if the Robertson problem -%# is formulated in explicit form (cf. odepkg_testsuite_implrober.m). -%# function mass = odepkg_testsuite_implrobermass (t, y, varargin) -%# mass = [1, 0, 0; 0, 1, 0; 0, 0, 0]; - -%# Returns the REFERENCE values for the implicit ROBERTSON problem -function y = odepkg_testsuite_implroberref () - y(1) = 0.20833401497012e-07; - y(2) = 0.83333607703347e-13; - y(3) = 0.99999997916650e+00; - -%!demo -%! vsolver = {@odebdi}; -%! for vcnt=1:length (vsolver) -%! virob{vcnt,1} = odepkg_testsuite_implrober (vsolver{vcnt}, 1e-7); -%! end -%! virob - -%# Local Variables: *** -%# mode: octave *** -%# End: *** +%# Copyright (C) 2006, Thomas Treichl <treichl@users.sourceforge.net> +%# OdePkg - Package for solving ordinary differential equations with octave +%# +%# 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 2 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, write to the Free Software +%# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + +%# -*- texinfo -*- +%# @deftypefn {Function File} {[@var{solution}] =} odepkg_testsuite_implrober (@var{@@solver}, @var{reltol}) +%# +%# 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). +%# +%# Run examples with the command +%# @example +%# demo odepkg_testsuite_implrober +%# @end example +%# @end deftypefn +%# +%# @seealso{odepkg} + +function vret = odepkg_testsuite_implrober (vhandle, vrtol) + + %# Check number and types of all input arguments + if (nargin ~= 2) + help ('odepkg_testsuite_implrober'); + error ('OdePkg:InvalidInputArgument', ... + 'Number of input arguments must be exactly two'); + elseif (~isa (vhandle, 'function_handle') || ~isscalar (vrtol)) + print_usage; + end + + vret{1} = vhandle; %# The name for the solver that is used + vret{2} = vrtol; %# The value for the realtive tolerance + vret{3} = vret{2} * 1e-2; %# The value for the absolute tolerance + vret{4} = vret{2}; %# The value for the first time step + %# Write a debug message on the screen, because this testsuite function + %# may be called more than once from a loop over all present solvers + fprintf (1, ['Testsuite implicit ROBERTSON, testing solver %7s with relative', ... + ' tolerance %2.0e\n'], func2str (vret{1}), vrtol); fflush (1); + + %# Setting the integration algorithm options + vstart = 0; %# The point of time when solving is started + vstop = 1e11; %# The point of time when solving is stoped + [vinity, vinityd] = odepkg_testsuite_implroberinit; %# The initial values + + vopt = odeset ('Refine', 0, 'RelTol', vret{2}, 'AbsTol', vret{3}, ... + 'InitialStep', vret{4}, 'Stats', 'on', 'NormControl', 'off', ... + 'Jacobian', @odepkg_testsuite_implroberjac, 'MaxStep', vstop-vstart); + %# 'OutputFcn', @odeplot); + + %# Calculate the algorithm, start timer and do solving + tic; vsol = feval (vhandle, @odepkg_testsuite_implroberfun, ... + [vstart, vstop], vinity, vinityd', vopt); + vret{12} = toc; %# The value for the elapsed time + vref = odepkg_testsuite_implroberref; %# Get the reference solution vector + if (max (size (vsol.y(end,:))) == max (size (vref))), vlst = vsol.y(end,:); + elseif (max (size (vsol.y(:,end))) == max (size (vref))), vlst = vsol.y(:,end); + end + vret{5} = odepkg_testsuite_calcmescd (vlst, vref, vret{3}, vret{2}); + vret{6} = odepkg_testsuite_calcscd (vlst, vref, vret{3}, vret{2}); + %# vret{7} = vsol.stats.success + vsol.stats.failed; %# The value for all evals + %# vret{8} = vsol.stats.success; %# The value for success evals + %# vret{9} = vsol.stats.fevals; %# The value for fun calls + %# vret{10} = vsol.stats.partial; %# The value for partial derivations + %# vret{11} = vsol.stats.ludecom; %# The value for LU decompositions + vret{7} = vsol.stats.nsteps + vsol.stats.nfailed; %# The value for all evals + vret{8} = vsol.stats.nsteps; %# The value for success evals + vret{9} = vsol.stats.nfevals; %# The value for fun calls + vret{10} = vsol.stats.nsolves; %# The value for partial derivations + vret{11} = vsol.stats.ndecomps; %# The value for LU decompositions + +%#function odepkg_testsuite_implrober () +%# A = odeset ('RelTol', 1e-4, ... %# proprietary ode15i needs 1e-6 to be stable +%# 'AbsTol', [1e-6, 1e-10, 1e-6], ... +%# 'Jacobian', @odepkg_testsuite_implroberjac); +%# [y0, yd0] = odepkg_testsuite_implroberinit; +%# odebdi (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0', A) +%# [y0, yd0] = odepkg_testsuite_implroberinit; +%# odebdi (@odepkg_testsuite_implroberfun, [0, 1e11], y0, yd0') + +%# Return the results for the for the implicit ROBERTSON problem +function res = odepkg_testsuite_implroberfun (t, y, yd, varargin) + res(1,1) = -0.04 * y(1) + 1e4 * y(2) * y(3) - yd(1); + res(2,1) = 0.04 * y(1) - 1e4 * y(2) * y(3) - 3e7 * y(2)^2 - yd(2); + res(3,1) = y(1) + y(2) + y(3) - 1; + +%# Return the INITIAL values for the implicit ROBERTSON problem +function [y0, yd0] = odepkg_testsuite_implroberinit () + y0 = [1, 0, 0]; + yd0 = [-4e-2, 4e-2, 0]; + +%# Return the JACOBIAN matrix for the implicit ROBERTSON problem +function [dfdy, dfdyd] = odepkg_testsuite_implroberjac (t, y, yd, varargin) + dfdy(1,1) = -0.04; + dfdy(1,2) = 1e4 * y(3); + dfdy(1,3) = 1e4 * y(2); + dfdy(2,1) = 0.04; + dfdy(2,2) = -1e4 * y(3) - 6e7 * y(2); + dfdy(2,3) = -1e4 * y(2); + dfdy(3,1) = 1; + dfdy(3,2) = 1; + dfdy(3,3) = 1; + + dfdyd(1,1) = -1; + dfdyd(2,2) = -1; + dfdyd(3,3) = 0; + +%# For the implicit form of the Robertson problem a mass matrix is not +%# allowed. This mass matrix is only needed if the Robertson problem +%# is formulated in explicit form (cf. odepkg_testsuite_implrober.m). +%# function mass = odepkg_testsuite_implrobermass (t, y, varargin) +%# mass = [1, 0, 0; 0, 1, 0; 0, 0, 0]; + +%# Return the REFERENCE values for the implicit ROBERTSON problem +function y = odepkg_testsuite_implroberref () + y(1) = 0.20833401497012e-07; + y(2) = 0.83333607703347e-13; + y(3) = 0.99999997916650e+00; + +%!demo +%! vsolver = {@odebdi}; +%! for vcnt=1:length (vsolver) +%! virob{vcnt,1} = odepkg_testsuite_implrober (vsolver{vcnt}, 1e-7); +%! end +%! virob + +%# Local Variables: *** +%# mode: octave *** +%# End: ***