c----------------------------------------------------------------------- c these are the subroutines for the c c Plate c ODE of dimension 80 c c----------------------------------------------------------------------- subroutine prob(problm,neqn,tbegin,tend,ijac,mljac,mujac) IMPLICIT double precision (A-H,O-Z) character*(*) problm integer neqn,ijac,mljac,mujac double precision tbegin,tend PARAMETER(MX=8,MY=5,MACHS1=2,MACHS2=4) COMMON /TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, + OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT problm = 'plate' neqn = 80 tbegin = 0d0 tend = 7.0d0 ijac = 1 mljac = neqn mujac = neqn C --- SET DEFAULT VALUES NX=MX NY=MY NACHS1=MACHS1 NACHS2=MACHS2 NXM1=NX-1 NYm1=NY-1 NDEMI=NX*NY OMEGA=1000.0D+0 STIFFN=100.0D+0 WEIGHT=200.0D+0 DENOM=NX+1 DELX=2.0D+0/DENOM USH4=1.0D+0/(DELX**4) FAC=STIFFN*USH4 return end c----------------------------------------------------------------------- subroutine init(neqn,y) integer neqn double precision y(neqn) integer i do 10 i=1,neqn y(i) = 0d0 10 continue return end c----------------------------------------------------------------------- SUBROUTINE feval(N, X, Y, DF) IMPLICIT real*8 (A-H,O-Z) DIMENSION Y(N), DF(N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT c ################ print c write(6,*)' dans fplate, n,x=',n,x c############## C -------- LA BOUCLE ------- DO 1 I=1,NX DO 1 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- DF(K)=Y(K+NDEMI) C ------ POINT CENTRAL --- UC=16.D0*Y(K) IF(I.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-1) END IF IF(I.LT.NX)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+1) END IF IF(J.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-NX) END IF IF(J.LT.NY)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+NX) END IF IF(I.GT.1 .AND.J.GT.1 )UC=UC+2.D0*Y(K-NX-1) IF(I.LT.NX.AND.J.GT.1 )UC=UC+2.D0*Y(K-NX+1) IF(I.GT.1 .AND.J.LT.NY)UC=UC+2.D0*Y(K+NX-1) IF(I.LT.NX.AND.J.LT.NY)UC=UC+2.D0*Y(K+NX+1) IF(I.GT.2)UC=UC+Y(K-2) IF(I.LT.NXM1)UC=UC+Y(K+2) IF(J.GT.2)UC=UC+Y(K-2*NX) IF(J.LT.NYM1)UC=UC+Y(K+2*NX) IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN XI=I*DELX FORCE=EXP(-5.D0*(X-XI-2.D0)**2)+EXP(-5.D0*(X-XI-5.D0)**2) ELSE FORCE=0.D0 END IF DF(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT 1 CONTINUE RETURN END c----------------------------------------------------------------------- subroutine jeval(N,X,Y,DFY,MEBND) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT c ################ print c write(6,*)' dans jac , n,x=',n,x c############## C------ METTRE A ZERO ------- DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C -------- LA BOUCLE ------- DO 2 I=1,NX DO 2 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- DFY(K,K+NDEMI)=1.D0 C ------ POINT CENTRAL --- DFY(K+NDEMI,K)=-FAC*16.D0 IF(I.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-1)=FAC*8.D0 END IF IF(I.LT.NX)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+1)=FAC*8.D0 END IF IF(J.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-NX)=FAC*8.D0 END IF IF(J.LT.NY)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+NX)=FAC*8.D0 END IF IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=-FAC*2.D0 IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=-FAC*2.D0 IF(I.GT.2)DFY(K+NDEMI,K-2)=-FAC IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=-FAC IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=-FAC IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=-FAC DFY(K+NDEMI,K+NDEMI)= -OMEGA 2 CONTINUE RETURN END c----------------------------------------------------------------------- subroutine solut(neqn,true) integer neqn double precision true(neqn) c c true(1)=0.000490143813851336 true(2)=0.000980081485560611 true(3)=0.001462893811482190 true(4)=0.001915822464411935 true(5)=0.002285152533727002 true(6)=0.002461353376688549 true(7)=0.002254597413097122 true(8)=0.001438312591933600 true(9)=0.000849025149228402 true(10)=0.001697885005625757 true(11)=0.002535239886068847 true(12)=0.003323989552181772 true(13)=0.003977902193560667 true(14)=0.004320231736082990 true(15)=0.004025679955083897 true(16)=0.002643206356123840 true(17)=0.000980287627702671 true(18)=0.001960162971121222 true(19)=0.002925787622964379 true(20)=0.003831644928823870 true(21)=0.004570305067454005 true(22)=0.004922706753377098 true(23)=0.004509194826194244 true(24)=0.002876625183867201 true(25)=0.000849025149228402 true(26)= 0.001697885005625757 true(27)=0.002535239886068847 true(28)=0.003323989552181772 true(29)=0.003977902193560667 true(30)=0.004320231736082990 true(31)=0.004025679955083897 true(32)=0.002643206356123840 true(33)=0.000490143813851336 true(34)=0.000980081485560611 true(35)=0.001462893811482190 true(36)=0.001915822464411935 true(37)=0.002285152533727002 true(38)=0.002461353376688549 true(39)=0.002254597413097122 true(40)=0.001438312591933600 true(41)=- 0.001177590304545409 true(42)=-0.002409005827992214 true(43)=-0.003722140831656533 true(44)=-0.005078780056048207 true(45)=-0.006302661811097914 true(46)=-0.006973399942926759 true(47)=-0.006394575120415784 true(48)=-0.003960464551310118 true(49)=-0.002040148244040460 true(50)=-0.004174829877953482 true(51)=-0.006456510337516159 true(52)=-0.008832503276738242 true(53)=-0.011029624807177369 true(54)=-0.012352389570141255 true(55)=-0.011524177328690540 true(56)=-0.007253301886026949 true(57)=-0.002355180609090818 true(58)=-0.004818011655984428 true(59)=-0.007444281663313065 true(60)=-0.010157560112096413 true(61)=-0.012605323622195828 true(62)=-0.013946799885853517 true(63)=-0.012789150240831569 true(64)=-0.007920929102620235 true(65)=-0.002040148244040460 true(66)=-0.004174829877953482 true(67)=-0.006456510337516159 true(68)=-0.008832503276738242 true(69)=-0.011029624807177369 true(70)=-0.012352389570141255 true(71)=-0.011524177328690540 true(72)=-0.007253301886026949 true(73)=-0.001177590304545409 true(74)=-0.002409005827992214 true(75)=-0.003722140831656533 true(76)=-0.005078780056048207 true(77)=-0.006302661811097914 true(78)=-0.006973399942926759 true(79)=-0.006394575120415784 true(80)=-0.003960464551310118 return end