c----------------------------------------------------------------------- c c Problem POLLUTION c ODE of dimension 20 c c----------------------------------------------------------------------- subroutine feval(neqn,t,y,dy,ierr,rpar,ipar) integer neqn,ierr,ipar double precision t,y(neqn),dy(neqn),rpar double precision k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13, = k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25 double precision r(25) k1=.35d0 k2=.266d2 k3=.123d5 k4=.86d-3 k5=.82d-3 k6=.15d5 k7=.13d-3 k8=.24d5 k9=.165d5 k10=.9d4 k11=.22d-1 k12=.12d5 k13=.188d1 k14=.163d5 k15=.48d7 k16=.35d-3 k17=.175d-1 k18=.1d9 k19=.444d12 k20=.124d4 k21=.21d1 k22=.578d1 k23=.474d-1 k24=.178d4 k25=.312d1 r( 1) = k1 *y( 1) r( 2) = k2 *y( 2)*y(4) r( 3) = k3 *y( 5)*y(2) r( 4) = k4 *y( 7) r( 5) = k5 *y( 7) r( 6) = k6 *y( 7)*y(6) r( 7) = k7 *y( 9) r( 8) = k8 *y( 9)*y(6) r( 9) = k9 *y(11)*y(2) r(10) = k10*y(11)*y(1) r(11) = k11*y(13) r(12) = k12*y(10)*y(2) r(13) = k13*y(14) r(14) = k14*y( 1)*y(6) r(15) = k15*y( 3) r(16) = k16*y( 4) r(17) = k17*y( 4) r(18) = k18*y(16) r(19) = k19*y(16) r(20) = k20*y(17)*y(6) r(21) = k21*y(19) r(22) = k22*y(19) r(23) = k23*y( 1)*y(4) r(24) = k24*y(19)*y(1) r(25) = k25*y(20) dy(1) = -r(1)-r(10)-r(14)-r(23)-r(24)+r(2)+r(3)+r(9)+r(11)+ = r(12)+r(22)+r(25) dy(2) = -r(2)-r(3)-r(9)-r(12)+r(1)+r(21) dy(3) = -r(15)+r(1)+r(17)+r(19)+r(22) dy(4) = -r(2)-r(16)-r(17)-r(23)+r(15) dy(5) = -r(3)+r(4)+r(4)+r(6)+r(7)+r(13)+r(20) dy(6) = -r(6)-r(8)-r(14)-r(20)+r(3)+r(18)+r(18) dy(7) = -r(4)-r(5)-r(6)+r(13) dy(8) = r(4)+r(5)+r(6)+r(7) dy(9) = -r(7)-r(8) dy(10) = -r(12)+r(7)+r(9) dy(11) = -r(9)-r(10)+r(8)+r(11) dy(12) = r(9) dy(13) = -r(11)+r(10) dy(14) = -r(13)+r(12) dy(15) = r(14) dy(16) = -r(18)-r(19)+r(16) dy(17) = -r(20) dy(18) = r(20) dy(19) = -r(21)-r(22)-r(24)+r(23)+r(25) dy(20) = -r(25)+r(24) return end c----------------------------------------------------------------------- subroutine init(neqn,y) integer neqn double precision y(neqn) y(1) = 0d0 y(2) = 0.2d0 y(3) = 0d0 y(4) = 0.04d0 y(5) = 0d0 y(6) = 0d0 y(7) = 0.1d0 y(8) = 0.3d0 y(9) = 0.01d0 y(10) = 0d0 y(11) = 0d0 y(12) = 0d0 y(13) = 0d0 y(14) = 0d0 y(15) = 0d0 y(16) = 0d0 y(17) = 0.007d0 y(18) = 0d0 y(19) = 0d0 y(20) = 0d0 return end c----------------------------------------------------------------------- subroutine jeval(neqn,t,y,jac,ldim,ierr,rpar,ipar) integer neqn,ldim,ierr,ipar double precision t,y(neqn),jac(ldim,neqn),rpar double precision k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13, = k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25 integer i,j do 20 j=1,neqn do 10 i=1,neqn jac(i,j)=0d0 10 continue 20 continue k1=.35d0 k2=.266d2 k3=.123d5 k4=.86d-3 k5=.82d-3 k6=.15d5 k7=.13d-3 k8=.24d5 k9=.165d5 k10=.9d4 k11=.22d-1 k12=.12d5 k13=.188d1 k14=.163d5 k15=.48d7 k16=.35d-3 k17=.175d-1 k18=.1d9 k19=.444d12 k20=.124d4 k21=.21d1 k22=.578d1 k23=.474d-1 k24=.178d4 k25=.312d1 jac(1,1) = -k1-k10*y(11)-k14*y(6)-k23*y(4)-k24*y(19) jac(1,11) = -k10*y(1)+k9*y(2) jac(1,6) = -k14*y(1) jac(1,4) = -k23*y(1)+k2*y(2) jac(1,19) = -k24*y(1)+k22 jac(1,2) = k2*y(4)+k9*y(11)+k3*y(5)+k12*y(10) jac(1,13) = k11 jac(1,20) = k25 jac(1,5) = k3*y(2) jac(1,10) = k12*y(2) jac(2,4) = -k2*y(2) jac(2,5) = -k3*y(2) jac(2,11) = -k9*y(2) jac(2,10) = -k12*y(2) jac(2,19) = k21 jac(2,1) = k1 jac(2,2) = -k2*y(4)-k3*y(5)-k9*y(11)-k12*y(10) jac(3,1) = k1 jac(3,4) = k17 jac(3,16) = k19 jac(3,19) = k22 jac(3,3) = -k15 jac(4,4) = -k2*y(2)-k16-k17-k23*y(1) jac(4,2) = -k2*y(4) jac(4,1) = -k23*y(4) jac(4,3) = k15 jac(5,5) = -k3*y(2) jac(5,2) = -k3*y(5) jac(5,7) = 2d0*k4+k6*y(6) jac(5,6) = k6*y(7)+k20*y(17) jac(5,9) = k7 jac(5,14) = k13 jac(5,17) = k20*y(6) jac(6,6) = -k6*y(7)-k8*y(9)-k14*y(1)-k20*y(17) jac(6,7) = -k6*y(6) jac(6,9) = -k8*y(6) jac(6,1) = -k14*y(6) jac(6,17) = -k20*y(6) jac(6,2) = k3*y(5) jac(6,5) = k3*y(2) jac(6,16) = 2d0*k18 jac(7,7) = -k4-k5-k6*y(6) jac(7,6) = -k6*y(7) jac(7,14) = k13 jac(8,7) = k4+k5+k6*y(6) jac(8,6) = k6*y(7) jac(8,9) = k7 jac(9,9) = -k7-k8*y(6) jac(9,6) = -k8*y(9) jac(10,10) = -k12*y(2) jac(10,2) = -k12*y(10)+k9*y(11) jac(10,9) = k7 jac(10,11) = k9*y(2) jac(11,11) = -k9*y(2)-k10*y(1) jac(11,2) = -k9*y(11) jac(11,1) = -k10*y(11) jac(11,9) = k8*y(6) jac(11,6) = k8*y(9) jac(11,13) = k11 jac(12,11) = k9*y(2) jac(12,2) = k9*y(11) jac(13,13) = -k11 jac(13,11) = k10*y(1) jac(13,1) = k10*y(11) jac(14,14) = -k13 jac(14,10) = k12*y(2) jac(14,2) = k12*y(10) jac(15,1) = k14*y(6) jac(15,6) = k14*y(1) jac(16,16) = -k18-k19 jac(16,4) = k16 jac(17,17) = -k20*y(6) jac(17,6) = -k20*y(17) jac(18,17) = k20*y(6) jac(18,6) = k20*y(17) jac(19,19) = -k21-k22-k24*y(1) jac(19,1) = -k24*y(19)+k23*y(4) jac(19,4) = k23*y(1) jac(19,20) = k25 jac(20,20) = -k25 jac(20,1) = k24*y(19) jac(20,19) = k24*y(1) return end c----------------------------------------------------------------------- subroutine prob(problm,neqn,ndisc,t,ijac,mljac,mujac) character*(*) problm integer neqn,ijac,mljac,mujac,ndisc double precision t(0:*) problm = 'pollut' neqn = 20 ndisc = 0 t(0) = 0d0 t(1) = 60d0 ijac = 1 mljac = neqn mujac = neqn return end c----------------------------------------------------------------------- subroutine solut(neqn,y) integer neqn double precision y(neqn) y( 1)= .5646255480022769d-1 y( 2)= .1342484130422339d+0 y( 3)= .4139734331099427d-8 y( 4)= .5523140207484359d-2 y( 5)= .2018977262302196d-6 y( 6)= .1464541863493966d-6 y( 7)= .7784249118997964d-1 y( 8)= .3245075353396018d+0 y( 9)= .7494013383880406d-2 y(10)= .1622293157301561d-7 y(11)= .1135863833257075d-7 y(12)= .2230505975721359d-2 y(13)= .2087162882798630d-3 y(14)= .1396921016840158d-4 y(15)= .8964884856898295d-2 y(16)= .4352846369330103d-17 y(17)= .6899219696263405d-2 y(18)= .1007803037365946d-3 y(19)= .1772146513969984d-5 y(20)= .5682943292316392d-4 return end