From octave-sources-request at bevo dot che dot wisc dot edu  Mon Jul 27 15:04:48 1998
Subject: Another test problem for diffeq
From: "Billinghurst, David (RTD)" <David dot Billinghurst at riotinto dot com dot au>
To: "'octave-sources at bevo dot che dot wisc dot edu'" 	 <octave-sources@bevo.che.wisc.edu>
Date: Mon, 27 Jul 1998 02:22:53 -0000

This message is in MIME format. Since your mail reader does not understand
this format, some or all of this message may not be legible.

------ =_NextPart_000_01BDB999.C909ACD2
Content-Type: text/plain

Here is another test problem for diffeq.  
The test case is problem CHEMAZKO from CWI Test Set for IVP Solvers
Ref: http://www.cwi.nl/cwi/projects/IVPtestset.shtml

+++++++++++++++++++++++++++++++++++++++++
(Mr) David Billinghurst
Comalco Research and Technical Support
PO Box 316, Thomastown, Vic, Australia, 3074
Phone:	+61 3 9469 0642
FAX:	+61 3 9462 2700
Email:	David dot Billinghurst at riotinto dot com

 <<dassl_chemakzo.m>>  <<lsode_chemakzo.m>> 

------ =_NextPart_000_01BDB999.C909ACD2
Content-Type: application/octet-stream;
	name="dassl_chemakzo.m"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
	filename="dassl_chemakzo.m"

# dassl_chemakzo.m=0A=
#=0A=
# problem CHEMAZKO  from CWI Test Set for IVP Solvers=0A=
# Ref: http://www.cwi.nl/cwi/projects/IVPtestset.shtml=0A=
#=0A=
# ODE of dimension 6=0A=
#=0A=
# Author: David Billinghurst (David dot Billinghurst at riotinto dot com dot au)=0A=
#         Comalco Research and Technology=0A=
#         Melbourne, Australia=0A=
#         22 July 1998=0A=
#=0A=
dassl_options("absolute tolerance", 1e-12);=0A=
dassl_options("relative tolerance", 1e-12);=0A=
tol =3D 1e-10;=0A=
=0A=
# Number of output points (including initial point)=0A=
# Note:  I could not achieve the required accuracy if N=3D2=0A=
N =3D 181;=0A=
Tout =3D 180;=0A=
=0A=
# The system of equations as dy/dt =3D f(y)=0A=
function dy =3D f(y)=0A=
  # Parameters=0A=
  k1   =3D 18.7;=0A=
  k2   =3D 0.58;=0A=
  k3   =3D 0.09;=0A=
  k4   =3D 0.42;=0A=
  K    =3D 34.4;=0A=
  klA  =3D 3.3;=0A=
  pCO2 =3D 0.9;=0A=
  H    =3D 737;=0A=
=0A=
  r1  =3D k1*(y(1)^4)*sqrt(y(2)); =0A=
  r2  =3D k2*y(3)*y(4);=0A=
  r3  =3D (k2/K)*y(1)*y(5);=0A=
  r4  =3D k3*y(1)*(y(4)^2);=0A=
  r5  =3D k4*(y(6)^2)*sqrt(y(2));=0A=
  Fin =3D klA*(pCO2/H-y(2));=0A=
=0A=
  dy =3D zeros(6,1);=0A=
  dy(1) =3D -2*r1 + r2 - r3 - r4;=0A=
  dy(2) =3D -0.5*r1  - r4 - 0.5*r5 + Fin;=0A=
  dy(3) =3D  r1 - r2 + r3;=0A=
  dy(4) =3D  - r2 + r3 - 2*r4;=0A=
  dy(5) =3D  r2 - r3 + r5;=0A=
  dy(6) =3D  - r5;=0A=
endfunction=0A=
=0A=
# The system of equations as g(y,ydot,t)=3D0=0A=
function res =3D g( y, ydot, t)=0A=
  res =3D ydot - f(y);=0A=
endfunction=0A=
=0A=
# Initial conditions=0A=
y0 =3D [ 0.437; 0.00123; 0; 0; 0; 0.367 ];=0A=
ydot0 =3D f(y0);=0A=
t =3D [0:Tout/(N-1):Tout];=0A=
=0A=
# Solution y at t =3D 180 seconds=0A=
soln =3D  [ 0.1161602274780192,=0A=
          0.001119418166040848,=0A=
          0.1621261719785814,=0A=
          0.003396981299297459,=0A=
          0.1646185108335055,=0A=
          0.1989533275954281 ];=0A=
=0A=
[y,ydot] =3D dassl("g",y0,ydot0,t);=0A=
y180 =3D y(N,:)';=0A=
=0A=
# Report "ans =3D 1" if final error is small (less than tol)=0A=
all(abs(y180-soln)<tol)=0A=
=0A=

------ =_NextPart_000_01BDB999.C909ACD2
Content-Type: application/octet-stream;
	name="lsode_chemakzo.m"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
	filename="lsode_chemakzo.m"

# lsode_chemakzo.m=0A=
#=0A=
# problem CHEMAZKO from CWI Test Set for IVP Solvers=0A=
# Ref: http://www.cwi.nl/cwi/projects/IVPtestset.shtml=0A=
#=0A=
# ODE of dimension 6=0A=
#=0A=
# Author: David Billinghurst (David dot Billinghurst at riotinto dot com dot au)=0A=
#         Comalco Research and Technology=0A=
#         Melbourne, Australia=0A=
#         22 July 1998=0A=
#=0A=
#=0A=
lsode_options("absolute tolerance", 1e-12);=0A=
lsode_options("relative tolerance", 1e-12);=0A=
tol =3D 1e-10;=0A=
=0A=
# Number of output points (including initial point)=0A=
# Note:  I could not achieve the required accuracy if N=3D2=0A=
N =3D 181;=0A=
Tout =3D 180;=0A=
=0A=
# The system of equations as dy/dt =3D f(y)=0A=
function dy =3D f(y)=0A=
  # Parameters=0A=
  k1   =3D 18.7;=0A=
  k2   =3D 0.58;=0A=
  k3   =3D 0.09;=0A=
  k4   =3D 0.42;=0A=
  K    =3D 34.4;=0A=
  klA  =3D 3.3;=0A=
  pCO2 =3D 0.9;=0A=
  H    =3D 737;=0A=
=0A=
  r1  =3D k1*(y(1)^4)*sqrt(y(2)); =0A=
  r2  =3D k2*y(3)*y(4);=0A=
  r3  =3D (k2/K)*y(1)*y(5);=0A=
  r4  =3D k3*y(1)*(y(4)^2);=0A=
  r5  =3D k4*(y(6)^2)*sqrt(y(2));=0A=
  Fin =3D klA*(pCO2/H-y(2));=0A=
=0A=
  dy =3D zeros(6,1);=0A=
  dy(1) =3D -2*r1 + r2 - r3 - r4;=0A=
  dy(2) =3D -0.5*r1  - r4 - 0.5*r5 + Fin;=0A=
  dy(3) =3D  r1 - r2 + r3;=0A=
  dy(4) =3D  -r2 + r3 - 2*r4;=0A=
  dy(5) =3D  r2 - r3 + r5;=0A=
  dy(6) =3D  -r5;=0A=
endfunction=0A=
=0A=
# Initial conditions=0A=
y0 =3D [ 0.437; 0.00123; 0; 0; 0; 0.367 ];=0A=
t =3D [0:Tout/(N-1):Tout];=0A=
=0A=
# Solution y at t =3D 180 seconds=0A=
soln =3D  [ 0.1161602274780192,=0A=
          0.001119418166040848,=0A=
          0.1621261719785814,=0A=
          0.003396981299297459,=0A=
          0.1646185108335055,=0A=
          0.1989533275954281 ];=0A=
=0A=
y =3D lsode("f",y0,t);=0A=
y180 =3D y(N,:)';=0A=
=0A=
# Report "ans =3D 1" if final error is small (less than tol)=0A=
all(abs(y180-soln)<tol)=0A=

------ =_NextPart_000_01BDB999.C909ACD2--