home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Monster Media 1994 #1
/
monster.zip
/
monster
/
PROG_C
/
C_LSODE.ZIP
/
TESTE.DOC
< prev
Wrap
Text File
|
1993-12-29
|
48KB
|
1,614 lines
teste.doc - test.doc for lsodE
this file contains five demos and their outputs. the file origin for
the first two files is /netlib/odepack/demo.Z at research.att.com and
for imsl.for is page 647-648, math/library, imsl manual, version 1.1,
jan, 1989. they are either trimmed or modified for lsode package.
contents
1. demoe.c
2. demoe.for
3. output of demoe.c
4. output of demoe.for
5. ivpag.for
6. imsl.for
7. imsl.c
8. outputs of 5, 6, and 7
9. outputs from example.for and example.c
(example.c is in lsode.h)
(example.for is in lsode.c)
**************************************************************************
**************************************************************************
********************************* demo.c *********************************
**************************************************************************
**************************************************************************
/* -------------------------------------------------------------------- */
/* demoe.c was translated from the file "demo" that was downloaded */
/* from netlib at research.att.com to test clsode. */
/* -------------------------------------------------------------------- */
/* demonstration program for the lsode package. */
/* this is the version of 13 august, 1981. */
/* */
/* this version is in double precision. */
/* */
/* the package is used to solve two simple problems, */
/* one with a full jacobian, the other with a banded jacobian, */
/* with all 8 of the appropriate values of mf in each case. */
/* if the errors are too large, or other difficulty occurs, */
/* a warning message is printed. all output is on unit lout = 6. */
/* -------------------------------------------------------------------- */
#include <stdlib.h>
#include <math.h>
#include "lsode.h"
#define pd(a,b) pd[nr*(b-1)+a]
void fcn1(int[], double, double[], double[]);
void jac1(int[], double, double[], int, int, double[], int);
void fcn2(int[], double, double[], double[]);
void jac2(int[], double, double[], int, int, double[], int);
double edit2(double[], double);
int main()
{
int i, iopar, iopt, iout, istate, itask, itol, *iwork,
leniw, lenrw, liw, lrw, mband, mf, miter, miter1, meth,
ml, mu, neq[2], nerr, nfe, nfea, nje, nout, nqu, nst;
double atole[2], dtout, er, erm, ero, hu, rtol[2], *rwork, t, tout,
tout1, *y;
tout1 = 1.39283880203e0;
dtout = 2.214773875e0;
neq[1] = 2;
nerr = 0;
itol = 1;
rtol[1] = 0.0e0;
atole[1] = 1.0e-6;
lrw = 697;
liw = 45;
iwork = ivector(liw);
rwork = dvector(lrw);
y = dvector(neq[1]);
iopt = 0;
/* */
/* first problem ----------------------------------------------------- */
/* */
nout = 4;
printf("\f\n\n demonstration program for clsode package\n\n\n\n");
printf("\n problem 1.. van der pol oscillator..");
printf("\n xdotdot - 3*(1 - x**2)*xdot + x = 0,");
printf(" x(0) = 2, xdot(0) = 0");
printf("\n neq = %2d",neq[1]);
printf("\n itol = %3d rtol = %10.1e atol = %10.1e",
itol, rtol[1], atole[1]);
for(meth=1;meth<=2;meth++) {
for(miter1=1;miter1<=4;miter1++) {
miter = miter1 - 1;
mf = 10*meth + miter;
printf("\n\n\n\n mf = %3d\n\n\n",mf);
printf(" t\t\t x\t\t xdot");
printf(" nq h\n");
t = 0.0e0;
y[1] = 2.0e0;
y[2] = 0.0e0;
itask = 1;
istate = 1;
tout = tout1;
ero = 0.0e0;
for(iout=1;iout<=nout;iout++) {
lsode(fcn1,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
iopt,rwork,lrw,iwork,liw,jac1,mf,stdout);
hu = rwork[11];
nqu = iwork[14];
printf("\n %15.5e %16.5e %14.3e %5d %14.3e",
t, y[1], y[2], nqu, hu);
if (istate < 0) break;
iopar = iout - 2*(iout/2);
if (iopar == 0) {
er = fabs(y[1])/atole[1];
ero = max(ero,er);
if (er >= 1000.0e0) {
printf("\n\n warning.. error ");
printf("exceeds 1000 * tolerance\n\n");
nerr++;
}
}
tout += dtout;
}
if (istate < 0) nerr++;
nst = iwork[11];
nfe = iwork[12];
nje = iwork[13];
lenrw = iwork[17];
leniw = iwork[18];
nfea = nfe;
if (miter == 2) nfea = nfe - neq[1]*nje;
if (miter == 3) nfea = nfe - nje;
printf("\n\n final statistics for this run..");
printf("\n rwork size = %4d iwork size = %4d",lenrw,leniw);
printf("\n number of steps = %5d",nst);
printf("\n number of f-s = %5d",nfe);
printf("\n (excluding j-s) = %5d",nfea);
printf("\n number of j-s = %5d",nje);
printf("\n error overrun = %10.2e",ero);
}
}
/* */
/* second problem ---------------------------------------------------- */
/* */
neq[1] = 25;
y = (double *) realloc(y, (unsigned)(neq[1]+1)*sizeof(double));
ml = 5;
mu = 0;
iwork[1] = ml;
iwork[2] = mu;
mband = ml + mu + 1;
nout = 5;
printf("\f\n\n demonstration program for clsode package\n\n\n\n");
printf("\n problem 2.. ydot = a * y , where,");
printf(" a is a banded lower triangular matrix");
printf("\n derived from 2-d advection pde");
printf("\n neq = %3d ml = %2d mu = %2d",neq[1],ml,mu);
printf("\n itol = %3d rtol = %10.1e atol = %10.1e",
itol, rtol[1], atole[1]);
for(meth=1;meth<=2;meth++) {
for(miter1=1;miter1<=6;miter1++) {
miter = miter1 - 1;
if (miter == 1 || miter == 2) goto skip_step;
mf = 10*meth + miter;
printf("\n\n\n\n mf = %3d\n\n\n",mf);
printf(" t\t\t max.err nq h\n");
t = 0.0e0;
for(i=2;i<=neq[1];i++) y[i] = 0.0e0;
y[1] = 1.0e0;
itask = 1;
istate = 1;
tout = 0.01e0;
ero = 0.0e0;
for(iout=1;iout<=nout;iout++) {
lsode(fcn2,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
iopt,rwork,lrw,iwork,liw,jac2,mf,stdout);
erm = edit2(y,t);
hu = rwork[11];
nqu = iwork[14];
printf("\n %15.5e %14.3e %5d %14.3e",t,erm,nqu,hu);
if (istate < 0) break;
er = erm/atole[1];
ero = max(ero,er);
if (er > 1000.0e0) {
printf("\n\n warning.. error ");
printf("exceeds 1000 * tolerance\n\n");
nerr++;
}
tout *= 10.0e0;
}
if (istate < 0) nerr++;
nst = iwork[11];
nfe = iwork[12];
nje = iwork[13];
lenrw = iwork[17];
leniw = iwork[18];
nfea = nfe;
if (miter == 5) nfea = nfe - mband*nje;
if (miter == 3) nfea = nfe - nje;
printf("\n\n final statistics for this run..");
printf("\n rwork size = %4d iwork size = %4d",lenrw,leniw);
printf("\n number of steps = %5d",nst);
printf("\n number of f-s = %5d",nfe);
printf("\n (excluding j-s) = %5d",nfea);
printf("\n number of j-s = %5d",nje);
printf("\n error overrun = %10.2e",ero);
skip_step: continue;
}
}
printf("\n\n\n\n number of errors encountered = %3d",nerr);
free(iwork);
free(rwork);
free(y);
return 1;
}
/* */
/* derivative and jacobi matrix routine ------------------------------ */
/* */
void fcn1(int neq[], double t, double y[], double ydot[])
{
ydot[1] = y[2];
ydot[2] = 3.0e0*(1.0e0 - y[1]*y[1])*y[2] - y[1];
return;
}
void jac1(int neq[], double t, double y[],
int ml, int mu, double pd[], int nr)
{
pd(1,1) = 0.0e0;
pd(1,2) = 1.0e0;
pd(2,1) = -6.0e0*y[1]*y[2] - 1.0e0;
pd(2,2) = 3.0e0*(1.0e0 - y[1]*y[1]);
return;
}
void fcn2(int neq[], double t, double y[], double ydot[])
{
int i, j, k, ng= 5;
double alph1 = 1.0e0, alph2 = 1.0e0, d;
for(j = 1;j<=ng;j++) {
for(i=1;i<=ng;i++) {
k = i + (j - 1)*ng;
d = -2.0e0*y[k];
if (i != 1) d = d + y[k-1]*alph1;
if (j != 1) d = d + y[k-ng]*alph2;
ydot[k] = d;
}
}
return;
}
void jac2(int neq[], double t, double y[],
int ml, int mu, double pd[], int nr)
{
int j, mband, mu1, mu2, ng=5;
double alph1 = 1.0e0, alph2 = 1.0e0;
mband = ml + mu + 1;
mu1 = mu + 1;
mu2 = mu + 2;
for(j=1;j<=neq[1];j++) {
pd(mu1,j) = -2.0e0;
pd(mu2,j) = alph1;
pd(mband,j) = alph2;
}
for(j=ng;j<=neq[1];j+=ng) pd(mu2,j) = 0.0e0;
return;
}
double edit2 (double y[], double t)
{
int i, j, k, ng=5;
double alph1=1.0e0, alph2=1.0e0;
double erm, a1, a2, er, ex, yt;
erm = 0.0e0;
if (t == 0.0e0) return erm;
ex = 0.0e0;
if (t <= 30.0e0) ex = exp(-2.0e0*t);
a2 = 1.0e0;
for(j=1;j<=ng;j++) {
a1 = 1.0e0;
for(i=1;i<=ng;i++) {
k = i + (j - 1)*ng;
yt = pow(t,(i+j-2))*ex*a1*a2;
er = fabs(y[k]-yt);
erm = max(erm,er);
a1 *= alph1/(double)i;
}
a2 *= alph2/(double)j;
}
return erm;
}
/* ---------------------------- end of demoe.c ------------------------ */
**************************************************************************
**************************************************************************
******************************* demoe.for ********************************
**************************************************************************
**************************************************************************
c-----------------------------------------------------------------------
c demonstration program for the lsode package.
c this is the version of 13 august, 1981.
c
c this version is in double precision.
c
c for computer systems requiring a program card, the following (with
c the c in column 1 removed) may be used..
c program lsdem(lsout,tape6=lsout)
c
c the package is used to solve two simple problems,
c one with a full jacobian, the other with a banded jacobian,
c with all 8 of the appropriate values of mf in each case.
c if the errors are too large, or other difficulty occurs,
c a warning message is printed. all output is on unit lout = 6.
c-----------------------------------------------------------------------
external f1, jac1, f2, jac2
integer i, iopar, iopt, iout, istate, itask, itol, iwork,
1 leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter, miter1,
2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t,
1 tout, tout1, y
dimension y(25), rwork(697), iwork(45)
data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/
c
nerr = 0
itol = 1
rtol = 0.0d0
atol = 1.0d-6
lrw = 697
liw = 45
iopt = 0
c
c first problem
c
neq = 2
nout = 4
write (lout,110) neq,itol,rtol,atol
110 format(1h1/1x,40h demonstration program for lsode package////
1 1x,39h problem 1.. van der pol oscillator../
2 1x,38h xdotdot - 3*(1 - x**2)*xdot + x = 0, ,
3 24h x(0) = 2, xdot(0) = 0/
4 1x,6h neq =,i2/
5 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,e10.1//)
c
do 195 meth = 1,2
do 190 miter1 = 1,4
miter = miter1 - 1
mf = 10*meth + miter
write (lout,120) mf
120 format(////1x,5h mf =,i3///
1 6x,1ht,15x,1hx,15x,4hxdot,7x,2hnq,6x,1hh/)
t = 0.0d0
y(1) = 2.0d0
y(2) = 0.0d0
itask = 1
istate = 1
tout = tout1
ero = 0.0d0
do 170 iout = 1,nout
call lsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
1 iopt,rwork,lrw,iwork,liw,jac1,mf)
hu = rwork(11)
nqu = iwork(14)
write (lout,140) t,y(1),y(2),nqu,hu
140 format(1x,e15.5,e16.5,e14.3,i5,e14.3)
if (istate .lt. 0) go to 175
iopar = iout - 2*(iout/2)
if (iopar .ne. 0) go to 170
er = dabs(y(1))/atol
ero = dmax1(ero,er)
if (er .lt. 1000.0d0) go to 170
write (lout,150)
150 format(//1x,41h warning.. error exceeds 1000 * tolerance//)
nerr = nerr + 1
170 tout = tout + dtout
175 continue
if (istate .lt. 0) nerr = nerr + 1
nst = iwork(11)
nfe = iwork(12)
nje = iwork(13)
lenrw = iwork(17)
leniw = iwork(18)
nfea = nfe
if (miter .eq. 2) nfea = nfe - neq*nje
if (miter .eq. 3) nfea = nfe - nje
write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
180 format(/1x,32h final statistics for this run../
1 1x,13h rwork size =,i4,15h iwork size =,i4/
2 1x,18h number of steps =,i5/
3 1x,18h number of f-s =,i5/
4 1x,18h (excluding j-s) =,i5/
5 1x,18h number of j-s =,i5/
6 1x,16h error overrun =,e10.2)
190 continue
195 continue
c
c second problem
c
neq = 25
ml = 5
mu = 0
iwork(1) = ml
iwork(2) = mu
mband = ml + mu + 1
nout = 5
write (lout,210) neq,ml,mu,itol,rtol,atol
210 format(1h1/1x,40h demonstration program for lsode package////
1 1x,33h problem 2.. ydot = a * y , where,
2 39h a is a banded lower triangular matrix/
3 1x,32h derived from 2-d advection pde/
4 1x,6h neq =,i3,7h ml =,i2,7h mu =,i2/
5 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,e10.1//)
do 295 meth = 1,2
do 290 miter1 = 1,6
miter = miter1 - 1
if (miter .eq. 1 .or. miter .eq. 2) go to 290
mf = 10*meth + miter
write (lout,220) mf
220 format(////1x,5h mf =,i3///
1 6x,1ht,13x,8hmax.err.,5x,2hnq,6x,1hh/)
t = 0.0d0
do 230 i = 2,neq
230 y(i) = 0.0d0
y(1) = 1.0d0
itask = 1
istate = 1
tout = 0.01d0
ero = 0.0d0
do 270 iout = 1,nout
call lsode(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
1 iopt,rwork,lrw,iwork,liw,jac2,mf)
call edit2(y,t,erm)
hu = rwork(11)
nqu = iwork(14)
write (lout,240) t,erm,nqu,hu
240 format(1x,e15.5,e14.3,i5,e14.3)
if (istate .lt. 0) go to 275
er = erm/atol
ero = dmax1(ero,er)
if (er .le. 1000.0d0) go to 270
write (lout,150)
nerr = nerr + 1
270 tout = tout*10.0d0
275 continue
if (istate .lt. 0) nerr = nerr + 1
nst = iwork(11)
nfe = iwork(12)
nje = iwork(13)
lenrw = iwork(17)
leniw = iwork(18)
nfea = nfe
if (miter .eq. 5) nfea = nfe - mband*nje
if (miter .eq. 3) nfea = nfe - nje
write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
290 continue
295 continue
write (lout,300) nerr
300 format(////1x,31h number of errors encountered =,i3)
stop
end
c
subroutine f1 (neq, t, y, ydot)
integer neq
double precision t, y, ydot
dimension y(2), ydot(2)
ydot(1) = y(2)
ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
return
end
subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
integer neq, ml, mu, nrowpd
double precision t, y, pd
dimension y(2), pd(nrowpd,2)
pd(1,1) = 0.0d0
pd(1,2) = 1.0d0
pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0
pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1))
return
end
c
subroutine f2 (neq, t, y, ydot)
integer neq, i, j, k, ng
double precision t, y, ydot, alph1, alph2, d
dimension y(1), ydot(1)
data alph1/1.0d0/, alph2/1.0d0/, ng/5/
do 10 j = 1,ng
do 10 i = 1,ng
k = i + (j - 1)*ng
d = -2.0d0*y(k)
if (i .ne. 1) d = d + y(k-1)*alph1
if (j .ne. 1) d = d + y(k-ng)*alph2
10 ydot(k) = d
return
end
c
subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
double precision t, y, pd, alph1, alph2
dimension y(1), pd(nrowpd,1)
data alph1/1.0d0/, alph2/1.0d0/, ng/5/
mband = ml + mu + 1
mu1 = mu + 1
mu2 = mu + 2
do 10 j = 1,neq
pd(mu1,j) = -2.0d0
pd(mu2,j) = alph1
10 pd(mband,j) = alph2
do 20 j = ng,neq,ng
20 pd(mu2,j) = 0.0d0
return
end
c
subroutine edit2 (y, t, erm)
integer i, j, k, ng
double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt
dimension y(25)
data alph1/1.0d0/, alph2/1.0d0/, ng/5/
erm = 0.0d0
if (t .eq. 0.0d0) return
ex = 0.0d0
if (t .le. 30.0d0) ex = dexp(-2.0d0*t)
a2 = 1.0d0
do 60 j = 1,ng
a1 = 1.0d0
do 50 i = 1,ng
k = i + (j - 1)*ng
yt = t**(i+j-2)*ex*a1*a2
er = dabs(y(k)-yt)
erm = dmax1(erm,er)
a1 = a1*alph1/dfloat(i)
50 continue
a2 = a2*alph2/dfloat(j)
60 continue
return
end
c ------------------------------ end of demo.for ---------------------
**************************************************************************
**************************************************************************
*************************** output of demoe.c ****************************
**************************************************************************
**************************************************************************
demonstration program for clsode package
problem 1.. van der pol oscillator..
xdotdot - 3*(1 - x**2)*xdot + x = 0, x(0) = 2, xdot(0) = 0
neq = 2
itol = 1 rtol = 0.0e+00 atol = 1.0e-06
mf = 10
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 3 1.229e-01
3.60761e+00 -7.79864e-05 -3.169e+00 5 2.171e-02
5.82239e+00 -1.68009e+00 2.911e-01 3 4.753e-02
8.03716e+00 1.16694e-04 3.169e+00 5 2.342e-02
final statistics for this run..
rwork size = 52 iwork size = 20
number of steps = 297
number of f-s = 352
(excluding j-s) = 352
number of j-s = 0
error overrun = 1.17e+02
mf = 11
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 1.205e-01
3.60761e+00 -1.77321e-05 -3.169e+00 5 1.872e-02
5.82239e+00 -1.68010e+00 2.911e-01 6 9.633e-02
8.03716e+00 2.58940e-05 3.169e+00 5 1.899e-02
final statistics for this run..
rwork size = 58 iwork size = 22
number of steps = 203
number of f-s = 281
(excluding j-s) = 281
number of j-s = 29
error overrun = 2.59e+01
mf = 12
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 1.205e-01
3.60761e+00 -1.77321e-05 -3.169e+00 5 1.872e-02
5.82239e+00 -1.68010e+00 2.911e-01 6 9.633e-02
8.03716e+00 2.58939e-05 3.169e+00 5 1.899e-02
final statistics for this run..
rwork size = 58 iwork size = 22
number of steps = 203
number of f-s = 339
(excluding j-s) = 281
number of j-s = 29
error overrun = 2.59e+01
mf = 13
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 7.386e-02
3.60761e+00 3.44007e-05 -3.169e+00 6 2.605e-02
5.82239e+00 -1.68011e+00 2.911e-01 4 1.328e-01
8.03716e+00 -5.91397e-05 3.169e+00 5 2.045e-02
final statistics for this run..
rwork size = 56 iwork size = 20
number of steps = 198
number of f-s = 315
(excluding j-s) = 289
number of j-s = 26
error overrun = 5.91e+01
mf = 20
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 5.490e-02
3.60761e+00 -5.65788e-05 -3.169e+00 5 1.429e-02
5.82239e+00 -1.68010e+00 2.911e-01 4 5.829e-02
8.03716e+00 1.03869e-04 3.169e+00 5 1.488e-02
final statistics for this run..
rwork size = 38 iwork size = 20
number of steps = 289
number of f-s = 321
(excluding j-s) = 321
number of j-s = 0
error overrun = 1.04e+02
mf = 21
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 6.762e-02
3.60761e+00 -4.89765e-05 -3.169e+00 5 1.413e-02
5.82239e+00 -1.68010e+00 2.911e-01 5 1.256e-01
8.03716e+00 9.68674e-05 3.169e+00 5 1.420e-02
final statistics for this run..
rwork size = 44 iwork size = 22
number of steps = 262
number of f-s = 345
(excluding j-s) = 345
number of j-s = 30
error overrun = 9.69e+01
mf = 22
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 6.762e-02
3.60761e+00 -4.89765e-05 -3.169e+00 5 1.413e-02
5.82239e+00 -1.68010e+00 2.911e-01 5 1.256e-01
8.03716e+00 9.68674e-05 3.169e+00 5 1.420e-02
final statistics for this run..
rwork size = 44 iwork size = 22
number of steps = 262
number of f-s = 405
(excluding j-s) = 345
number of j-s = 30
error overrun = 9.69e+01
mf = 23
t x xdot nq h
1.39284e+00 1.68010e+00 -2.911e-01 5 7.089e-02
3.60761e+00 -4.67048e-05 -3.169e+00 5 1.390e-02
5.82239e+00 -1.68010e+00 2.911e-01 3 7.192e-02
8.03716e+00 5.47001e-05 3.169e+00 5 1.538e-02
final statistics for this run..
rwork size = 42 iwork size = 20
number of steps = 271
number of f-s = 414
(excluding j-s) = 383
number of j-s = 31
error overrun = 5.47e+01
demonstration program for clsode package
problem 2.. ydot = a * y , where, a is a banded lower triangular matrix
derived from 2-d advection pde
neq = 25 ml = 5 mu = 0
itol = 1 rtol = 0.0e+00 atol = 1.0e-06
mf = 10
t max.err nq h
1.00000e-02 5.563e-07 2 7.660e-03
1.00000e-01 6.550e-06 3 2.494e-02
1.00000e+00 2.744e-06 4 5.198e-02
1.00000e+01 1.143e-06 3 1.169e-01
1.00000e+02 2.215e-06 2 2.615e-01
final statistics for this run..
rwork size = 420 iwork size = 20
number of steps = 524
number of f-s = 552
(excluding j-s) = 552
number of j-s = 0
error overrun = 6.55e+00
mf = 13
t max.err nq h
1.00000e-02 8.390e-07 2 9.493e-03
1.00000e-01 2.075e-06 3 2.497e-02
1.00000e+00 1.268e-04 3 1.677e-02
1.00000e+01 1.130e-05 3 3.853e-01
1.00000e+02 1.475e-06 2 1.486e+01
final statistics for this run..
rwork size = 447 iwork size = 20
number of steps = 129
number of f-s = 235
(excluding j-s) = 201
number of j-s = 34
error overrun = 1.27e+02
mf = 14
t max.err nq h
1.00000e-02 8.772e-07 2 9.648e-03
1.00000e-01 2.064e-06 3 2.497e-02
1.00000e+00 1.263e-06 5 9.348e-02
1.00000e+01 3.108e-07 6 4.423e-01
1.00000e+02 1.594e-08 2 2.909e+01
final statistics for this run..
rwork size = 697 iwork size = 45
number of steps = 92
number of f-s = 113
(excluding j-s) = 113
number of j-s = 18
error overrun = 2.06e+00
mf = 15
t max.err nq h
1.00000e-02 8.771e-07 2 9.648e-03
1.00000e-01 2.064e-06 3 2.497e-02
1.00000e+00 1.263e-06 5 9.348e-02
1.00000e+01 3.108e-07 6 4.423e-01
1.00000e+02 1.595e-08 2 2.909e+01
final statistics for this run..
rwork size = 697 iwork size = 45
number of steps = 92
number of f-s = 221
(excluding j-s) = 113
number of j-s = 18
error overrun = 2.06e+00
mf = 20
t max.err nq h
1.00000e-02 4.648e-07 2 4.827e-03
1.00000e-01 1.307e-06 3 1.477e-02
1.00000e+00 4.267e-06 5 6.347e-02
1.00000e+01 1.921e-06 4 3.514e-01
1.00000e+02 9.884e-08 1 4.526e-01
final statistics for this run..
rwork size = 245 iwork size = 20
number of steps = 330
number of f-s = 530
(excluding j-s) = 530
number of j-s = 0
error overrun = 4.27e+00
mf = 23
t max.err nq h
1.00000e-02 1.006e-06 2 5.983e-03
1.00000e-01 4.456e-07 3 1.459e-02
1.00000e+00 1.527e-06 5 7.378e-02
1.00000e+01 5.778e-07 4 3.235e-01
1.00000e+02 5.610e-12 1 4.747e+01
final statistics for this run..
rwork size = 272 iwork size = 20
number of steps = 189
number of f-s = 343
(excluding j-s) = 290
number of j-s = 53
error overrun = 1.53e+00
mf = 24
t max.err nq h
1.00000e-02 1.039e-06 2 6.081e-03
1.00000e-01 4.632e-07 3 1.459e-02
1.00000e+00 2.475e-06 5 6.658e-02
1.00000e+01 8.280e-07 5 3.912e-01
1.00000e+02 3.843e-10 1 1.084e+02
final statistics for this run..
rwork size = 522 iwork size = 45
number of steps = 118
number of f-s = 136
(excluding j-s) = 136
number of j-s = 18
error overrun = 2.47e+00
mf = 25
t max.err nq h
1.00000e-02 1.039e-06 2 6.081e-03
1.00000e-01 4.632e-07 3 1.459e-02
1.00000e+00 2.475e-06 5 6.658e-02
1.00000e+01 8.280e-07 5 3.912e-01
1.00000e+02 3.843e-10 1 1.084e+02
final statistics for this run..
rwork size = 522 iwork size = 45
number of steps = 118
number of f-s = 244
(excluding j-s) = 136
number of j-s = 18
error overrun = 2.47e+00
number of errors encountered = 0
**************************************************************************
**************************************************************************
************************** output of demoe.for ***************************
**************************************************************************
**************************************************************************
demonstration program for lsode package
problem 1.. van der pol oscillator..
xdotdot - 3*(1 - x**2)*xdot + x = 0, x(0) = 2, xdot(0) = 0
neq = 2
itol = 1 rtol = .0E+00 atol = .1E-05
mf = 10
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 3 .123E+00
.36076E+01 -.77986E-04 -.317E+01 5 .217E-01
.58224E+01 -.16801E+01 .291E+00 3 .475E-01
.80372E+01 .11669E-03 .317E+01 5 .234E-01
final statistics for this run..
rwork size = 52 iwork size = 20
number of steps = 297
number of f-s = 352
(excluding j-s) = 352
number of j-s = 0
error overrun = .12E+03
mf = 11
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .121E+00
.36076E+01 -.17732E-04 -.317E+01 5 .187E-01
.58224E+01 -.16801E+01 .291E+00 6 .963E-01
.80372E+01 .25894E-04 .317E+01 5 .190E-01
final statistics for this run..
rwork size = 58 iwork size = 22
number of steps = 203
number of f-s = 281
(excluding j-s) = 281
number of j-s = 29
error overrun = .26E+02
mf = 12
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .121E+00
.36076E+01 -.17732E-04 -.317E+01 5 .187E-01
.58224E+01 -.16801E+01 .291E+00 6 .963E-01
.80372E+01 .25894E-04 .317E+01 5 .190E-01
final statistics for this run..
rwork size = 58 iwork size = 22
number of steps = 203
number of f-s = 339
(excluding j-s) = 281
number of j-s = 29
error overrun = .26E+02
mf = 13
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .739E-01
.36076E+01 .34401E-04 -.317E+01 6 .260E-01
.58224E+01 -.16801E+01 .291E+00 4 .133E+00
.80372E+01 -.59140E-04 .317E+01 5 .205E-01
final statistics for this run..
rwork size = 56 iwork size = 20
number of steps = 198
number of f-s = 315
(excluding j-s) = 289
number of j-s = 26
error overrun = .59E+02
mf = 20
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .549E-01
.36076E+01 -.56579E-04 -.317E+01 5 .143E-01
.58224E+01 -.16801E+01 .291E+00 4 .583E-01
.80372E+01 .10387E-03 .317E+01 5 .149E-01
final statistics for this run..
rwork size = 38 iwork size = 20
number of steps = 289
number of f-s = 321
(excluding j-s) = 321
number of j-s = 0
error overrun = .10E+03
mf = 21
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .676E-01
.36076E+01 -.48977E-04 -.317E+01 5 .141E-01
.58224E+01 -.16801E+01 .291E+00 5 .126E+00
.80372E+01 .96867E-04 .317E+01 5 .142E-01
final statistics for this run..
rwork size = 44 iwork size = 22
number of steps = 262
number of f-s = 345
(excluding j-s) = 345
number of j-s = 30
error overrun = .97E+02
mf = 22
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .676E-01
.36076E+01 -.48977E-04 -.317E+01 5 .141E-01
.58224E+01 -.16801E+01 .291E+00 5 .126E+00
.80372E+01 .96867E-04 .317E+01 5 .142E-01
final statistics for this run..
rwork size = 44 iwork size = 22
number of steps = 262
number of f-s = 405
(excluding j-s) = 345
number of j-s = 30
error overrun = .97E+02
mf = 23
t x xdot nq h
.13928E+01 .16801E+01 -.291E+00 5 .709E-01
.36076E+01 -.46705E-04 -.317E+01 5 .139E-01
.58224E+01 -.16801E+01 .291E+00 3 .719E-01
.80372E+01 .54700E-04 .317E+01 5 .154E-01
final statistics for this run..
rwork size = 42 iwork size = 20
number of steps = 271
number of f-s = 414
(excluding j-s) = 383
number of j-s = 31
error overrun = .55E+02
demonstration program for lsode package
problem 2.. ydot = a * y , where a is a banded lower triangular matrix
derived from 2-d advection pde
neq = 25 ml = 5 mu = 0
itol = 1 rtol = .0E+00 atol = .1E-05
mf = 10
t max.err. nq h
.10000E-01 .556E-06 2 .766E-02
.10000E+00 .655E-05 3 .249E-01
.10000E+01 .274E-05 4 .520E-01
.10000E+02 .114E-05 3 .117E+00
.10000E+03 .221E-05 2 .262E+00
final statistics for this run..
rwork size = 420 iwork size = 20
number of steps = 524
number of f-s = 552
(excluding j-s) = 552
number of j-s = 0
error overrun = .65E+01
mf = 13
t max.err. nq h
.10000E-01 .839E-06 2 .949E-02
.10000E+00 .208E-05 3 .250E-01
.10000E+01 .127E-03 3 .168E-01
.10000E+02 .113E-04 3 .385E+00
.10000E+03 .147E-05 2 .149E+02
final statistics for this run..
rwork size = 447 iwork size = 20
number of steps = 129
number of f-s = 235
(excluding j-s) = 201
number of j-s = 34
error overrun = .13E+03
mf = 14
t max.err. nq h
.10000E-01 .877E-06 2 .965E-02
.10000E+00 .206E-05 3 .250E-01
.10000E+01 .126E-05 5 .935E-01
.10000E+02 .311E-06 6 .442E+00
.10000E+03 .159E-07 2 .291E+02
final statistics for this run..
rwork size = 697 iwork size = 45
number of steps = 92
number of f-s = 113
(excluding j-s) = 113
number of j-s = 18
error overrun = .21E+01
mf = 15
t max.err. nq h
.10000E-01 .877E-06 2 .965E-02
.10000E+00 .206E-05 3 .250E-01
.10000E+01 .126E-05 5 .935E-01
.10000E+02 .311E-06 6 .442E+00
.10000E+03 .160E-07 2 .291E+02
final statistics for this run..
rwork size = 697 iwork size = 45
number of steps = 92
number of f-s = 221
(excluding j-s) = 113
number of j-s = 18
error overrun = .21E+01
mf = 20
t max.err. nq h
.10000E-01 .465E-06 2 .483E-02
.10000E+00 .131E-05 3 .148E-01
.10000E+01 .427E-05 5 .635E-01
.10000E+02 .192E-05 4 .351E+00
.10000E+03 .988E-07 1 .453E+00
final statistics for this run..
rwork size = 245 iwork size = 20
number of steps = 330
number of f-s = 530
(excluding j-s) = 530
number of j-s = 0
error overrun = .43E+01
mf = 23
t max.err. nq h
.10000E-01 .101E-05 2 .598E-02
.10000E+00 .446E-06 3 .146E-01
.10000E+01 .153E-05 5 .738E-01
.10000E+02 .578E-06 4 .324E+00
.10000E+03 .561E-11 1 .475E+02
final statistics for this run..
rwork size = 272 iwork size = 20
number of steps = 189
number of f-s = 343
(excluding j-s) = 290
number of j-s = 53
error overrun = .15E+01
mf = 24
t max.err. nq h
.10000E-01 .104E-05 2 .608E-02
.10000E+00 .463E-06 3 .146E-01
.10000E+01 .247E-05 5 .666E-01
.10000E+02 .828E-06 5 .391E+00
.10000E+03 .384E-09 1 .108E+03
final statistics for this run..
rwork size = 522 iwork size = 45
number of steps = 118
number of f-s = 136
(excluding j-s) = 136
number of j-s = 18
error overrun = .25E+01
mf = 25
t max.err. nq h
.10000E-01 .104E-05 2 .608E-02
.10000E+00 .463E-06 3 .146E-01
.10000E+01 .247E-05 5 .666E-01
.10000E+02 .828E-06 5 .391E+00
.10000E+03 .384E-09 1 .108E+03
final statistics for this run..
rwork size = 522 iwork size = 45
number of steps = 118
number of f-s = 244
(excluding j-s) = 136
number of j-s = 18
error overrun = .25E+01
number of errors encountered = 0
Stop - Program terminated.
**************************************************************************
**************************************************************************
******************************* ivpag.for ********************************
**************************************************************************
**************************************************************************
C IVPAG.FOR
C FILE ORIGIN:
C IMSL USER'S MANUAL, IMSL MATH/LIBRARY, VERSION 1.1, JAN, 1989
C PAGE 647 - 648
C
INTEGER NEQ, NPARAM
PARAMETER (NEQ=3, NPARAM=50)
C
INTEGER IDO, IEND, IMETH, INORM, NOUT
REAL A(1,1), FCN, FCNJ, HINIT, PARAM(NPARAM), TOL, X,
& XEND, Y(NEQ)
EXTERNAL FCN, FCNJ, IVPAG, SSET, UMACH
C Initialize
HINIT = 1.0E-03
INORM = 2
IMETH = 1
CALL SSET (NPARAM, 0.0, PARAM, 1)
PARAM(1) = HINIT
PARAM(10) = INORM
PARAM(12) = IMETH
C
IDO = 1
X = 0.0
Y(1) = 0.0
Y(2) = 1.0
Y(3) = 1.0
TOL = 1.0E-6
C Write title
CALL UMACH (2, NOUT)
WRITE (NOUT,99998)
C Integrate ODE
DO 10 IEND=1, 10
XEND = IEND
CALL IVPAG (IDO, NEQ, FCN, FCNJ, A, X, XEND, TOL, PARAM, Y)
WRITE (NOUT,99999) X, Y
10 CONTINUE
C Finish up
IDO = 3
CALL IVPAG (IDO, NEQ, FCN, FCNJ, A, X, XEND, TOL, PARAM, Y)
C
99998 FORMAT (11X, 'X', 14X, 'Y(1)', 11X, 'Y(2)', 11X, 'Y(3)')
99999 FORMAT (4F15.5)
END
C
SUBROUTINE FCN (NEQ, X, Y, YPRIME)
INTEGER NEQ
DIMENSION Y(NEQ), YPRIME(NEQ)
C
YPRIME(1) = Y(2)*Y(3)
YPRIME(2) = -Y(1)*Y(3)
YPRIME(3) = -0.51D0*Y(1)*Y(2)
RETURN
END
C
SUBROUTINE FCNJ(NEQ, X, Y, DYPDY)
INTEGER NEQ
DIMENSION X, Y(NEQ), DYPDY(*)
C This subroutine is never called
RETURN
END
**************************************************************************
**************************************************************************
********************************* imsl.for *******************************
**************************************************************************
**************************************************************************
c
c imsl.for
c file origin:
c imsl user's manual, imsl math/library, version 1.1, jan, 1989
c page 647 - 648
c
double precision atol,rtol,rwork,t,tout,y
integer i,iopt,istate,itask,itol,iwork,liw,lrw,mf,neq
dimension iwork(23),rwork(68),y(3)
external fcn,jac
c
c connect unit 6 for output to the screen.
c
open (6,file='con')
c
c enable display of error messages on unit 6.
c
call xsetf(1)
call xsetun(6)
c
c initialization of lsode call list parameters.
c
neq = 3
y(1) = 0.0d0
y(2) = 1.0d0
y(3) = 1.0d0
c
t = 0.d0
tout = 0.0d0
c
itol = 1
rtol = 1.d-4
atol = 1.d-06
c
itask = 1
istate = 1
iopt = 1
lrw = 68
liw = 23
c
mf = 10
write(6,24)
c
c solution and output of results.
c
do 10, i = 0,10
call lsode(fcn,neq,y,t,tout,itol,rtol,atol,itask,istate,
* iopt,rwork,lrw,iwork,liw,fcnj,mf)
write(6,20) t,y
if (istate .lt. 0) then
write (6,60) istate
write (6,30) iwork(11)
write (6,40) iwork(12)
stop
endif
tout = tout + 1.0d0
10 continue
c
write (6,30) iwork(11)
write (6,40) iwork(12)
c
20 format (4f15.5)
24 format (/,6x,'time',14x,'y(1)',16x,'y(2)',16x,'y(3)'/)
25 format (3x,f9.5,6x,3(f12.7,6x))
30 format (/,1x,'number of steps = ',i4)
40 format (1x,'number of fcn evaluations =',i4)
60 format (/,1x,'error halt: istate = ',i2)
c
end
c
c
c fcn defines the system of odes.
c
subroutine fcn(neq, x, y, yprime)
c
double precision x, y, yprime
integer neq
dimension y(neq), yprime(neq)
c
yprime(1) = y(2)*y(3)
yprime(2) = -y(1)*y(3)
yprime(3) = -0.51d0*y(1)*y(2)
return
end
c
c
c jac computes the jacobian matrix.
c this subroutine is never called.
c
subroutine fcnj(neq,x,y,ml,mu,dypdy,nr)
c
double precision x, y, dypdy
integer ml,mu,neq,nr
dimension y(1), dypdy(nr,1)
c
return
end
**************************************************************************
**************************************************************************
********************************** imsl.c ********************************
**************************************************************************
**************************************************************************
/* imsl.c */
/* file origin: */
/* imsl user's manual, imsl math/library, version 1.1, jan, 1989 */
/* page 647 - 648 */
#include <math.h>
#include <stdlib.h>
#include "lsode.h"
void fcn(int[], double, double[], double[]);
int main()
{
double atole[1+1], rtol[1+1], *rwork,t,tout,*y;
int i,iopt,istate,itask,itol,*iwork,liw,lrw,mf,neq[1+1];
neq[1] = 3;
y = dvector(neq[1]);
y[1] = 0.0e0;
y[2] = 1.0e0;
y[3] = 1.0e0;
t = 0.0e0;
tout = 0.0e0;
itol = 1;
rtol[1] = 1.e-04;
atole[1] = 1.e-06;
itask = 1;
istate = 1;
iopt = 1;
lrw = 68;
liw = 23;
rwork = dvector(lrw);
iwork = ivector(liw);
mf = 10;
/* */
/* integration step - call lsode. ----------------------------------- */
/* */
for(i=0; i<=10; i++) {
lsode(fcn,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
iopt,rwork,lrw,iwork,liw,NULL,mf,stdout);
if (istate < 0) {
printf("\n error halt: istate = %3d",istate);
break;
}
printf("\n%15.5f %15.5f %15.5f %15.5f",t,y[1],y[2],y[3]);
tout += 1.0e0;
}
printf("\n\n number of steps = %d",iwork[11]);
printf("\n number of fcn evaluations = %d",iwork[12]);
free(iwork);
free(rwork);
free(y);
return 0;
}
void fcn(int neq[], double t, double y[], double yprime[])
{
yprime[1] = y[2]*y[3];
yprime[2] = -y[1]*y[3];
yprime[3] = -0.51e0*y[1]*y[2];
return;
}
**************************************************************************
**************************************************************************
**************** outputs of ivpag.for, imsl.for, and imsl.c **************
**************************************************************************
**************************************************************************
OUTPUT OF IVPAG.FOR (SINGLE PRECISION RESULT)
X Y(1) Y(2) Y(3)
1.00000 0.80220 0.59705 0.81963
2.00000 0.99537 -0.09615 0.70336
3.00000 0.64141 -0.76720 0.88892
4.00000 -0.26961 -0.96296 0.98129
5.00000 -0.91172 -0.41079 0.75899
6.00000 -0.95750 0.28841 0.72967
7.00000 -0.42877 0.90342 0.95197
8.00000 0.51093 0.85963 0.93106
9.00000 0.97567 0.21925 0.71730
10.00000 0.87790 -0.47885 0.77906
output of imsl.for (double precision result) n-s = 82 n-f = 89
time y(1) y(2) y(3)
1.00000 .80220 .59710 .81966
2.00000 .99528 -.09618 .70341
3.00000 .64130 -.76714 .88895
4.00000 -.26932 -.96319 .98143
5.00000 -.91187 -.41084 .75898
6.00000 -.95751 .28848 .72971
7.00000 -.42891 .90333 .95195
8.00000 .51058 .86030 .93137
9.00000 .97596 .21953 .71729
10.00000 .87823 -.47878 .77905
output of imsl.c (double precision result) n-s = 82 n-f = 89
time y[1] y[2] y[3]
1.00000 0.80220 0.59710 0.81966
2.00000 0.99528 -0.09618 0.70341
3.00000 0.64130 -0.76714 0.88895
4.00000 -0.26932 -0.96319 0.98143
5.00000 -0.91187 -0.41084 0.75898
6.00000 -0.95751 0.28848 0.72971
7.00000 -0.42891 0.90333 0.95195
8.00000 0.51058 0.86030 0.93137
9.00000 0.97596 0.21953 0.71729
10.00000 0.87823 -0.47878 0.77905