slsode
所属分类:数学计算
开发工具:Fortran
文件大小:78KB
下载次数:5
上传日期:2011-03-08 15:55:49
上 传 者:
ultrafast
说明: 可求解一般微分方程组的数学算法,高效率fortran求解器
(numerical algorithm)
文件列表:
slsode\Debug\DF60.PDB (33792, 2008-07-14)
slsode\Debug\slsode.obj (131834, 2008-07-14)
slsode\Debug\slsode.pdb (25600, 2008-07-14)
slsode\slsode.dsp (4153, 2008-07-14)
slsode\slsode.dsw (535, 2008-07-14)
slsode\slsode.for (170894, 2008-07-14)
slsode\slsode.ncb (41984, 2008-07-14)
slsode\slsode.opt (43520, 2008-07-14)
slsode\slsode.plg (600, 2008-07-14)
slsode\Debug (0, 2008-07-16)
slsode (0, 2008-07-16)
c-----------------------------------------------------------------------
c Demonstration program for the SLSODE package.
c This is the version of 12 June 2001.
c
c This version is in single precision.
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,
2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
real 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.39283880203e0/, dtout/2.214773875e0/
c
nerr = 0
itol = 1
rtol = 0.0e0
atol = 1.0e-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(/' Demonstration program for SLSODE package'///
1 ' Problem 1: Van der Pol oscillator:'/
2 ' xdotdot - 3*(1 - x**2)*xdot + x = 0, ',
3 ' x(0) = 2, xdot(0) = 0'/
4 ' neq =',i2/
5 ' itol =',i3,' rtol =',e10.1,' atol =',e10.1//)
c
do 195 meth = 1,2
do 190 miter = 0,3
mf = 10*meth + miter
write (lout,120) mf
120 format(///' Solution with mf =',i3//
1 5x,'t x xdot nq h'//)
t = 0.0e0
y(1) = 2.0e0
y(2) = 0.0e0
itask = 1
istate = 1
tout = tout1
ero = 0.0e0
do 170 iout = 1,nout
call slsode(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(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 = abs(y(1))/atol
ero = max(ero,er)
if (er .gt. 1000.0e0) then
write (lout,150)
150 format(//' Warning: error exceeds 1000 * tolerance'//)
nerr = nerr + 1
endif
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(//' Final statistics for this run:'/
1 ' rwork size =',i4,' iwork size =',i4/
2 ' number of steps =',i5/
3 ' number of f-s =',i5/
4 ' (excluding J-s) =',i5/
5 ' number of J-s =',i5/
6 ' 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(///70('-')///
1 ' Problem 2: ydot = A * y , where',
2 ' A is a banded lower triangular matrix'/
3 12x, 'derived from 2-D advection PDE'/
4 ' neq =',i3,' ml =',i2,' mu =',i2/
5 ' itol =',i3,' rtol =',e10.1,' atol =',e10.1//)
do 295 meth = 1,2
do 290 miter = 0,5
if (miter .eq. 1 .or. miter .eq. 2) go to 290
mf = 10*meth + miter
write (lout,220) mf
220 format(///' Solution with mf =',i3//
1 5x,'t max.err. nq h'//)
t = 0.0e0
do 230 i = 2,neq
230 y(i) = 0.0e0
y(1) = 1.0e0
itask = 1
istate = 1
tout = 0.01e0
ero = 0.0e0
do 270 iout = 1,nout
call slsode(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(e15.5,e14.3,i5,e14.3)
if (istate .lt. 0) go to 275
er = erm/atol
ero = max(ero,er)
if (er .gt. 1000.0e0) then
write (lout,150)
nerr = nerr + 1
endif
270 tout = tout*10.0e0
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(////' Number of errors encountered =',i3)
stop
end
subroutine f1 (neq, t, y, ydot)
integer neq
real t, y, ydot
dimension y(neq), ydot(neq)
ydot(1) = y(2)
ydot(2) = 3.0e0*(1.0e0 - y(1)*y(1))*y(2) - y(1)
return
end
subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
integer neq, ml, mu, nrowpd
real t, y, pd
dimension y(neq), pd(nrowpd,neq)
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
end
subroutine f2 (neq, t, y, ydot)
integer neq, i, j, k, ng
real t, y, ydot, alph1, alph2, d
dimension y(neq), ydot(neq)
data alph1/1.0e0/, alph2/1.0e0/, ng/5/
do 10 j = 1,ng
do 10 i = 1,ng
k = i + (j - 1)*ng
d = -2.0e0*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
subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
real t, y, pd, alph1, alph2
dimension y(neq), pd(nrowpd,neq)
data alph1/1.0e0/, alph2/1.0e0/, ng/5/
mband = ml + mu + 1
mu1 = mu + 1
mu2 = mu + 2
do 10 j = 1,neq
pd(mu1,j) = -2.0e0
pd(mu2,j) = alph1
10 pd(mband,j) = alph2
do 20 j = ng,neq,ng
20 pd(mu2,j) = 0.0e0
return
end
subroutine edit2 (y, t, erm)
integer i, j, k, ng
real y, t, erm, alph1, alph2, a1, a2, er, ex, yt
dimension y(25)
data alph1/1.0e0/, alph2/1.0e0/, ng/5/
erm = 0.0e0
if (t .eq. 0.0e0) return
ex = 0.0e0
if (t .le. 30.0e0) ex = exp(-2.0e0*t)
a2 = 1.0e0
do 60 j = 1,ng
a1 = 1.0e0
do 50 i = 1,ng
k = i + (j - 1)*ng
yt = t**(i+j-2)*ex*a1*a2
er = abs(y(k)-yt)
erm = max(erm,er)
a1 = a1*alph1/i
50 continue
a2 = a2*alph2/j
60 continue
return
end
Demonstration program for DLSODE 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 = 0.1E-05
Solution with mf = 10
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 3 0.123E+00
0.36076E+01 -0.77***6E-04 -0.317E+01 5 0.217E-01
0.58224E+01 -0.16801E+01 0.291E+00 3 0.475E-01
0.80372E+01 0.11669E-03 0.317E+01 5 0.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 = 0.12E+03
Solution with mf = 11
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00
0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01
0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01
0.80372E+01 0.25894E-04 0.317E+01 5 0.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 = 0.26E+02
Solution with mf = 12
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00
0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01
0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01
0.80372E+01 0.25894E-04 0.317E+01 5 0.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 = 0.26E+02
Solution with mf = 13
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.739E-01
0.36076E+01 0.34401E-04 -0.317E+01 6 0.260E-01
0.58224E+01 -0.16801E+01 0.291E+00 4 0.133E+00
0.80372E+01 -0.59053E-04 0.317E+01 5 0.205E-01
Final statistics for this run:
rwork size = 56 iwork size = 20
number of steps = 1***
number of f-s = 315
(excluding J-s) = 289
number of J-s = 26
error overrun = 0.59E+02
Solution with mf = 20
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.549E-01
0.36076E+01 -0.56579E-04 -0.317E+01 5 0.143E-01
0.58224E+01 -0.16801E+01 0.291E+00 4 0.583E-01
0.80372E+01 0.10387E-03 0.317E+01 5 0.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 = 0.10E+03
Solution with mf = 21
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01
0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01
0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00
0.80372E+01 0.96867E-04 0.317E+01 5 0.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 = 0.97E+02
Solution with mf = 22
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01
0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01
0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00
0.80372E+01 0.96867E-04 0.317E+01 5 0.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 = 0.97E+02
Solution with mf = 23
t x xdot nq h
0.13928E+01 0.16801E+01 -0.291E+00 5 0.709E-01
0.36076E+01 -0.46705E-04 -0.317E+01 5 0.139E-01
0.58224E+01 -0.16801E+01 0.291E+00 3 0.719E-01
0.80372E+01 0.54700E-04 0.317E+01 5 0.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 = 0.55E+02
----------------------------------------------------------------------
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 = 0.1E-05
Solution with mf = 10
t max.err. nq h
0.10000E-01 0.556E-06 2 0.766E-02
0.10000E+00 0.655E-05 3 0.249E-01
0.10000E+01 0.274E-05 4 0.520E-01
0.10000E+02 0.114E-05 3 0.117E+00
0.10000E+03 0.221E-05 2 0.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 = 0.65E+01
Solution with mf = 13
t max.err. nq h
0.10000E-01 0.839E-06 2 0.949E-02
0.10000E+00 0.208E-05 3 0.250E-01
0.10000E+01 0.127E-03 3 0.168E-01
0.10000E+02 0.113E-04 3 0.385E+00
0.10000E+03 0.145E-05 2 0.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 = 0.13E+03
Solution with mf = 14
t max.err. nq h
0.10000E-01 0.877E-06 2 0.965E-02
0.10000E+00 0.206E-05 3 0.250E-01
0.10000E+01 0.126E-05 5 0.935E-01
0.10000E+02 0.311E-06 6 0.442E+00
0.10000E+03 0.159E-07 2 0.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 = 0.21E+01
Solution with mf = 15
t max.err. nq h
0.10000E-01 0.877E-06 2 0.965E-02
0.10000E+00 0.206E-05 3 0.250E-01
0.10000E+01 0.126E-05 5 0.935E-01
0.10000E+02 0.311E-06 6 0.442E+00
0.10000E+03 0.160E-07 2 0.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 = 0.21E+01
Solution with mf = 20
t max.err. nq h
0.10000E-01 0.465E-06 2 0.483E-02
0.10000E+00 0.131E-05 3 0.148E-01
0.10000E+01 0.427E-05 5 0.635E-01
0.10000E+02 0.192E-05 4 0.351E+00
0.10000E+03 0.929E-07 1 0.455E+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 = 0.43E+01
Solution with mf = 23
t max.err. nq h
0.10000E-01 0.101E-05 2 0.5***E-02
0.10000E+00 0.446E-06 3 0.146E-01
0.10000E+01 0.153E-05 5 0.738E-01
0.10000E+02 0.578E-06 4 0.324E+00
0.10000E+03 0.908E-08 1 0.992E+02
Final statistics for this run:
rwork size = 272 iwork size = 20
number of steps = 180
number of f-s = 325
(excluding J-s) = 274
number of J-s = 51
error overrun = 0.15E+01
Solution with mf = 24
t max.err. nq h
0.10000E-01 0.104E-05 2 0.608E-02
0.10000E+00 0.463E-06 3 0.146E-01
0.10000E+01 0.247E-05 5 0.666E-01
0.10000E+02 0.828E-06 5 0.391E+00
0.10000E+03 0.384E-09 1 0.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 = 0.25E+01
Solution with mf = 25
t max.err. nq h
0.10000E-01 0.104E-05 2 0.608E-02
0.10000E+00 0.463E-06 3 0.146E-01
0.10000E+01 0.247E-05 5 0.666E-01
0.10000E+02 0.828E-06 5 0.391E+00
0.10000E+03 0.384E-09 1 0.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 = 0.25E+01
Number of errors encountered = 0
近期下载者:
相关文件:
收藏者: