C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT ANDREW PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * C LINK dr_andrew radau5m vecteurw decsol dc_decsol OR C LINK dr_andrew radau5m vecteurw lapack lapackc dc_lapack IMPLICIT DOUBLE PRECISION (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=20,LWORK=4*ND*ND+23*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) DIMENSION RPAR(41) CHARACTER*20 NOMFICY,NOMFICZ EXTERNAL FEX,JEX,SOLOUT,SOLOUTZ,MAS C --- DEFINITION OF THE RESULT FILES WRITE(6,*)'NAME OF THE RESULT FILE FOR THE DIFFERENTIAL COMPONENTS?' READ(*,*) NOMFICY WRITE(6,*)'THE RESULT FILE FOR Y IS : ',NOMFICY NFICY = 10 OPEN(NFICY,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICY) WRITE(6,*)'NAME OF THE RESULT FILE FOR THE ALGEBRAIC COMPONENTS?' READ(*,*) NOMFICZ WRITE(6,*)'THE RESULT FILE FOR Z IS : ',NOMFICZ NFICZ = 11 OPEN(NFICZ,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICZ) C --- PARAMETER IN THE DIFFERENTIAL EQUATION IPAR=1 C --- CONSTANT M1 RPAR(1) =.04325D0 C --- CONSTANT M2 RPAR(2) =.00365D0 C --- CONSTANT M3 RPAR(3) =.02373D0 C --- CONSTANT M4 RPAR(4) =.00706D0 C --- CONSTANT M5 RPAR(5) =.07050D0 C --- CONSTANT M6 RPAR(6) =.00706D0 C --- CONSTANT M7 RPAR(7) =.05498D0 C --- CONSTANT XB RPAR(8) =-0.03635D0 C --- CONSTANT YB RPAR(9) =.03273D0 C --- CONSTANT XC RPAR(10) =.014D0 C --- CONSTANT YC RPAR(11) =.072D0 C --- CONSTANT C0 RPAR(12) =4530.0D0 C --- CONSTANT I1 RPAR(13) =2.194D-6 C --- CONSTANT I2 RPAR(14) =4.410D-7 C --- CONSTANT I3 RPAR(15) =5.255D-6 C --- CONSTANT I4 RPAR(16) =5.667D-7 C --- CONSTANT I5 RPAR(17) =1.169D-5 C --- CONSTANT I6 RPAR(18) =5.667D-7 C --- CONSTANT I7 RPAR(19) =1.912D-5 C --- CONSTANT D RPAR(20) =0.028D0 C --- CONSTANT DA RPAR(21) =115.0D-4 C --- CONSTANT E RPAR(22) =2.0D-2 C --- CONSTANT EA RPAR(23) =1421.0D-5 C --- CONSTANT RR RPAR(24) =7.0D-3 C --- CONSTANT RA RPAR(25) =92.0D-5 C --- CONSTANT L0 RPAR(26) =7785.0D-5 C --- CONSTANT SS RPAR(27) =35.0D-3 C --- CONSTANT SA RPAR(28) =1874.0D-5 C --- CONSTANT SB RPAR(29) =1043.0D-5 C --- CONSTANT SC RPAR(30) =18.0D-3 C --- CONSTANT SD RPAR(31) =2.0D-2 C --- CONSTANT TA RPAR(32) =2308.0D-5 C --- CONSTANT TB RPAR(33) =916.0D-5 C --- CONSTANT U RPAR(34) =4.0D-2 C --- CONSTANT UA RPAR(35) =1228.0D-5 C --- CONSTANT UB RPAR(36) =449.0D-5 C --- CONSTANT ZF RPAR(37) =2.0D-2 C --- CONSTANT ZT RPAR(38) =4.0D-2 C --- CONSTANT FA RPAR(39) =1421.0D-5 C --- CONSTANT MOM RPAR(40) =33.0D-3 C --- CONSTANT PI RPAR(41) =4.0D0*DATAN(1.0D0) C --- DIMENSION OF THE SYSTEM N=20 C --- COMPUTE THE JACOBIAN NUMERICALY IJAC=0 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN IMPLICIT FORM IMAS=1 MLMAS=0 MUMAS=0 C --- OUTPUT ROUTINE IS USED DURING INTEGRATION IOUT=1 C --- INITIAL VALUES X=0.0D0 Y(1 )=-0.0617138900142764496358948458001D0 Y(2 )= 0D0 Y(3 )= 0.455279819163070380255912382449D0 Y(4 )= 0.222668390165885884674473185609D0 Y(5 )= 0.487364979543842550225598953530D0 Y(6 )= -0.222668390165885884674473185609D0 Y(7 )= 1.23054744454982119249735015568D0 Y(8 )= 0D0 Y(9 )= 0D0 Y(10)= 0D0 Y(11)= 0D0 Y(12)= 0D0 Y(13)= 0D0 Y(14)= 0D0 Y(15)= 98.5668703962410896057654982170D0 Y(16)= -6.12268834425566265503114393122D0 Y(17)= 0D0 Y(18)= 0D0 Y(19)= 0D0 Y(20)= 0D0 C --- ENDPOINT OF INTEGRATION XEND=3.0D-2 C --- REQUIRED TOLERANCE WRITE(*,*)'RTOL?' READ(*,*)RTOL ATOL=1.0D0*RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-7 C --- SET DEFAULT VALUES DO I=1,20 IWORK(I)=0 WORK(I)=0.D0 END DO IWORK(5)=14 IWORK(6)=6 C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5M(N,FEX,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JEX,IJAC,MLJAC,MUJAC, & MAS,IMAS,MLMAS,MUMAS, & SOLOUT,SOLOUTZ,IOUT,NFICY,NFICZ, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- CLOSE THE RESULT FILES CLOSE(NFICY) CLOSE(NFICZ) C --- PRINT FINAL SOLUTION WRITE (6,99) X,Y(1),Y(2) WRITE (6,99) X,Y(3),Y(4) WRITE (6,99) X,Y(5),Y(6) WRITE (6,99) X,Y(7),Y(8) WRITE (6,99) X,Y(9),Y(10) WRITE (6,99) X,Y(11),Y(12) WRITE (6,99) X,Y(13),Y(14) WRITE (6,98) X,Y(15),Y(16) WRITE (6,98) X,Y(17),Y(18) WRITE (6,98) X,Y(19),Y(20) 99 FORMAT(1X,'X =',F5.2,' Y =',2E22.14) 98 FORMAT(1X,'X =',F5.2,' Z =',2E22.14) C --- PRINT STATISTICS WRITE (6,90) RTOL 90 FORMAT(' RTOL=',D8.2) WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' FCN=',I5,' JAC=',I4,' STEP=',I4,' ACCPT=',I4, & ' REJCT=',I3,' DEC=',I4,' SOL=',I5) STOP END C C SUBROUTINE SOLOUT (NFICY,NFICZ,NR,XOLD,X,Y,CONT,LRC,N, & RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),CONT(LRC),RPAR(*) COMMON /INTERN/XOUTY,XOUTZ PI2=2*RPAR(41) IF (NR.EQ.1) THEN WRITE (NFICY,99) X,DMOD(Y(1),PI2),DMOD(Y(2),PI2),NR-1 WRITE (NFICY,99) X,DMOD(Y(3),PI2),DMOD(Y(4),PI2),NR-1 WRITE (NFICY,99) X,DMOD(Y(5),PI2),DMOD(Y(6),PI2),NR-1 WRITE (NFICY,99) X,DMOD(Y(7),PI2),Y(8),NR-1 WRITE (NFICY,99) X,Y(9),Y(10),NR-1 WRITE (NFICY,99) X,Y(11),Y(12),NR-1 WRITE (NFICY,99) X,Y(13),Y(14),NR-1 WRITE (NFICZ,99) X,Y(15),Y(16),NR-1 WRITE (NFICZ,99) X,Y(17),Y(18),NR-1 WRITE (NFICZ,99) X,Y(19),Y(20),NR-1 XOUTY=3.0D-4 ELSE 10 CONTINUE IF (X.GE.XOUTY) THEN C --- CONTINUOUS OUTPUT FOR RADAU5M WRITE (NFICY,99) XOUTY,DMOD(CONTR5(1,XOUTY,CONT,LRC),PI2), & DMOD(CONTR5(2,XOUTY,CONT,LRC),PI2),NR-1 WRITE (NFICY,99) XOUTY,DMOD(CONTR5(3,XOUTY,CONT,LRC),PI2), & DMOD(CONTR5(4,XOUTY,CONT,LRC),PI2),NR-1 WRITE (NFICY,99) XOUTY,DMOD(CONTR5(5,XOUTY,CONT,LRC),PI2), & DMOD(CONTR5(6,XOUTY,CONT,LRC),PI2),NR-1 WRITE (NFICY,99) XOUTY,DMOD(CONTR5(7,XOUTY,CONT,LRC),PI2), & CONTR5(8,XOUTY,CONT,LRC),NR-1 WRITE (NFICY,99) XOUTY,CONTR5(9,XOUTY,CONT,LRC), & CONTR5(10,XOUTY,CONT,LRC),NR-1 WRITE (NFICY,99) XOUTY,CONTR5(11,XOUTY,CONT,LRC), & CONTR5(12,XOUTY,CONT,LRC),NR-1 WRITE (NFICY,99) XOUTY,CONTR5(13,XOUTY,CONT,LRC), & CONTR5(14,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,99) XOUTY,CONTR5(15,XOUTY,CONT,LRC), & CONTR5(16,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,99) XOUTY,CONTR5(17,XOUTY,CONT,LRC), & CONTR5(18,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,99) XOUTY,CONTR5(19,XOUTY,CONT,LRC), & CONTR5(20,XOUTY,CONT,LRC),NR-1 XOUTY=XOUTY+3.0D-4 GOTO 10 END IF END IF 99 FORMAT(1X,F8.6,2E22.14,I4) RETURN END C C SUBROUTINE SOLOUTZ (NFICY,NFICZ,NR,XOLD,X,Y,N,VALG,LRV, & RPAR,IPAR,IRTRN,ICAS,R1,R2,HTOT,H2,H3,LAST) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "VECTWY" AND "VECTWZ" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),RPAR(*) DIMENSION VWT(9),WY(6),VALG(LRV) COMMON/INTERN/XOUTY,XOUTZ EXTERNAL VECTWZ,VECTWY PI2 = 2*RPAR(41) IF (NR.EQ.1) THEN WRITE(NFICZ,99) X,Y(15),Y(16),NR-1 WRITE(NFICZ,99) X,Y(17),Y(18),NR-1 WRITE(NFICZ,99) X,Y(19),Y(20),NR-1 WRITE(NFICY,99) X,DMOD(Y(1),PI2),DMOD(Y(2),PI2),NR-1 WRITE(NFICY,99) X,DMOD(Y(3),PI2),DMOD(Y(4),PI2),NR-1 WRITE(NFICY,99) X,DMOD(Y(5),PI2),DMOD(Y(6),PI2),NR-1 WRITE(NFICY,99) X,DMOD(Y(7),PI2),Y(8),NR-1 WRITE(NFICY,99) X,Y(9),Y(10),NR-1 WRITE(NFICY,99) X,Y(11),Y(12),NR-1 WRITE(NFICY,99) X,Y(13),Y(14),NR-1 XOUTZ=3.0D-4 XOUTY=3.0D-4 ELSE IF (LAST) THEN XPAR=X ELSE XPAR=XOLD ENDIF 51 CONTINUE IF (XPAR.GE.XOUTZ) THEN T=(XOUTZ+HTOT-X)/HTOT CALL VECTWZ(ICAS,VWT,R1,R2,T) W1=VWT(1) W2=VWT(2) W3=VWT(3) W4=VWT(4) W5=VWT(5) W6=VWT(6) W7=VWT(7) W8=VWT(8) W9=VWT(9) YN15=W1*VALG(15)+W2*VALG(75)+W3*VALG(135)+ + W4*VALG(35)+W5*VALG(95)+W6*VALG(155)+ + W7*VALG(55)+W8*VALG(115)+W9*VALG(175) YN16=W1*VALG(16)+W2*VALG(76)+W3*VALG(136)+ + W4*VALG(36)+W5*VALG(96)+W6*VALG(156)+ + W7*VALG(56)+W8*VALG(116)+W9*VALG(176) YN17=W1*VALG(17)+W2*VALG(77)+W3*VALG(137)+ + W4*VALG(37)+W5*VALG(97)+W6*VALG(157)+ + W7*VALG(57)+W8*VALG(117)+W9*VALG(177) YN18=W1*VALG(18)+W2*VALG(78)+W3*VALG(138)+ + W4*VALG(38)+W5*VALG(98)+W6*VALG(158)+ + W7*VALG(58)+W8*VALG(118)+W9*VALG(178) YN19=W1*VALG(19)+W2*VALG(79)+W3*VALG(139)+ + W4*VALG(39)+W5*VALG(99)+W6*VALG(159)+ + W7*VALG(59)+W8*VALG(119)+W9*VALG(179) YN20=W1*VALG(20)+W2*VALG(80)+W3*VALG(140)+ + W4*VALG(40)+W5*VALG(100)+W6*VALG(160)+ + W7*VALG(60)+W8*VALG(120)+W9*VALG(180) WRITE(NFICZ,99) XOUTZ,YN15,YN16,NR-1 WRITE(NFICZ,99) XOUTZ,YN17,YN18,NR-1 WRITE(NFICZ,99) XOUTZ,YN19,YN20,NR-1 XOUTZ=XOUTZ+3.0D-4 GOTO 51 ENDIF 52 CONTINUE IF (X.GE.XOUTY) THEN IF ((XOLD-H2).GE.XOUTY) THEN HTOTY=HTOT-H3 T=(XOUTY-XOLD+HTOTY)/HTOTY R=(HTOT-H3-H2)/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(1)+WY2*VALG(61)+WY3*VALG(121) YN1=YN1+WY4*VALG(21)+WY5*VALG(81)+WY6*VALG(141) YN2=WY1*VALG(2)+WY2*VALG(62)+WY3*VALG(122) YN2=YN2+WY4*VALG(22)+WY5*VALG(82)+WY6*VALG(142) YN3=WY1*VALG(3)+WY2*VALG(63)+WY3*VALG(123) YN3=YN3+WY4*VALG(23)+WY5*VALG(83)+WY6*VALG(143) YN4=WY1*VALG(4)+WY2*VALG(64)+WY3*VALG(124) YN4=YN4+WY4*VALG(24)+WY5*VALG(84)+WY6*VALG(144) YN5=WY1*VALG(5)+WY2*VALG(65)+WY3*VALG(125) YN5=YN5+WY4*VALG(25)+WY5*VALG(125)+WY6*VALG(145) YN6=WY1*VALG(6)+WY2*VALG(66)+WY3*VALG(126) YN6=YN6+WY4*VALG(26)+WY5*VALG(86)+WY6*VALG(146) YN7=WY1*VALG(7)+WY2*VALG(67)+WY3*VALG(127) YN7=YN7+WY4*VALG(27)+WY5*VALG(87)+WY6*VALG(147) YN8=WY1*VALG(8)+WY2*VALG(68)+WY3*VALG(128) YN8=YN8+WY4*VALG(28)+WY5*VALG(88)+WY6*VALG(148) YN9=WY1*VALG(9)+WY2*VALG(69)+WY3*VALG(129) YN9=YN9+WY4*VALG(29)+WY5*VALG(89)+WY6*VALG(149) YN10=WY1*VALG(10)+WY2*VALG(70)+WY3*VALG(130) YN10=YN10+WY4*VALG(30)+WY5*VALG(90)+WY6*VALG(150) YN11=WY1*VALG(11)+WY2*VALG(71)+WY3*VALG(131) YN11=YN11+WY4*VALG(31)+WY5*VALG(91)+WY6*VALG(151) YN12=WY1*VALG(12)+WY2*VALG(72)+WY3*VALG(132) YN12=YN12+WY4*VALG(32)+WY5*VALG(92)+WY6*VALG(152) YN13=WY1*VALG(13)+WY2*VALG(73)+WY3*VALG(133) YN13=YN13+WY4*VALG(33)+WY5*VALG(93)+WY6*VALG(153) YN14=WY1*VALG(14)+WY2*VALG(74)+WY3*VALG(134) YN14=YN14+WY4*VALG(34)+WY5*VALG(94)+WY6*VALG(154) ELSE HTOTY=H2+H3 T=(XOUTY-X+HTOTY)/HTOTY R=H2/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(21)+WY2*VALG(81)+WY3*VALG(141) YN1=YN1+WY4*VALG(41)+WY5*VALG(101)+WY6*VALG(161) YN2=WY1*VALG(22)+WY2*VALG(82)+WY3*VALG(142) YN2=YN2+WY4*VALG(42)+WY5*VALG(102)+WY6*VALG(162) YN3=WY1*VALG(23)+WY2*VALG(83)+WY3*VALG(143) YN3=YN3+WY4*VALG(43)+WY5*VALG(103)+WY6*VALG(163) YN4=WY1*VALG(24)+WY2*VALG(84)+WY3*VALG(144) YN4=YN4+WY4*VALG(44)+WY5*VALG(104)+WY6*VALG(164) YN5=WY1*VALG(25)+WY2*VALG(85)+WY3*VALG(145) YN5=YN5+WY4*VALG(45)+WY5*VALG(105)+WY6*VALG(165) YN6=WY1*VALG(26)+WY2*VALG(86)+WY3*VALG(146) YN6=YN6+WY4*VALG(46)+WY5*VALG(106)+WY6*VALG(166) YN7=WY1*VALG(27)+WY2*VALG(87)+WY3*VALG(147) YN7=YN7+WY4*VALG(47)+WY5*VALG(107)+WY6*VALG(167) YN8=WY1*VALG(28)+WY2*VALG(88)+WY3*VALG(148) YN8=YN8+WY4*VALG(48)+WY5*VALG(108)+WY6*VALG(168) YN9=WY1*VALG(29)+WY2*VALG(89)+WY3*VALG(149) YN9=YN9+WY4*VALG(49)+WY5*VALG(109)+WY6*VALG(169) YN10=WY1*VALG(30)+WY2*VALG(90)+WY3*VALG(150) YN10=YN10+WY4*VALG(50)+WY5*VALG(110)+WY6*VALG(170) YN11=WY1*VALG(31)+WY2*VALG(91)+WY3*VALG(151) YN11=YN11+WY4*VALG(51)+WY5*VALG(111)+WY6*VALG(171) YN12=WY1*VALG(32)+WY2*VALG(92)+WY3*VALG(152) YN12=YN12+WY4*VALG(52)+WY5*VALG(112)+WY6*VALG(172) YN13=WY1*VALG(33)+WY2*VALG(93)+WY3*VALG(153) YN13=YN13+WY4*VALG(53)+WY5*VALG(113)+WY6*VALG(173) YN14=WY1*VALG(34)+WY2*VALG(94)+WY3*VALG(154) YN14=YN14+WY4*VALG(54)+WY5*VALG(114)+WY6*VALG(174) ENDIF WRITE(NFICY,99) XOUTY,DMOD(YN1,PI2), & DMOD(YN2,PI2),NR-1 WRITE(NFICY,99) XOUTY,DMOD(YN3,PI2), & DMOD(YN4,PI2),NR-1 WRITE(NFICY,99) XOUTY,DMOD(YN5,PI2), & DMOD(YN6,PI2),NR-1 WRITE(NFICY,99) XOUTY,DMOD(YN7,PI2),YN8,NR-1 WRITE(NFICY,99) XOUTY,YN9,YN10,NR-1 WRITE(NFICY,99) XOUTY,YN11,YN12,NR-1 WRITE(NFICY,99) XOUTY,YN13,YN14,NR-1 XOUTY=XOUTY+3.0D-4 GOTO 52 ENDIF ENDIF RETURN 99 FORMAT(1X,F8.6,2E22.14,I4) END C C SUBROUTINE INVMAT(N,A,AINV,E,IP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(N,*),AINV(N,*),E(*) INTEGER IP(*) EXTERNAL DEC,SOL C --- DECOMPOSITION LU OF THE MATRIX A CALL DEC(N,N,A,IP,IER) IF (IER.NE.0) GOTO 10 C --- COMPUTAION OF THE INVERSE MATRIX DO I=1,N DO J=1,N E(J)=0.0D0 ENDDO E(I)=1.0D0 CALL SOL(N,N,A,E,IP) DO K=1,N AINV(K,I)=E(K) ENDDO ENDDO RETURN 10 WRITE(*,*)'SINGULAR MATRIX IN INVMAT',IER STOP END C C SUBROUTINE FEX(N,X,Y,DY,RPAR,IPAR) C --- RIGHT-HAND SIDE OF PENDULUM PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION M1,M2,M3,M4,M5,M6,M7,XB,YB,XC,YC,C0, + I1,I2,I3,I4,I5,I6,I7,D,DA,E,EA,RR,RA,L0, + SS,SA,SB,SC,SD,TA,TB,U,UA,UB,ZF,ZT,FA,MOM, + SIBE,SITH,SIGA,SIPH,SIDE,SIOM,SIEP, + COBE,COTH,COGA,COPH,CODE,COOM,COEP, + SIBETH,SIPHDE,SIOMEP,COBETH,COPHDE,COOMEP, + BEP,THP,PHP,DEP,OMP,EPP, + M,F,GP,XD,YD,LANG,FORCE,FX,FY,MINV,EI INTEGER IP,NM DIMENSION Y(N),DY(N),RPAR(*) DIMENSION M(7,7),F(7),GP(6,7),MINV(7,7),EI(7),IP(7) M1 = RPAR(1) M2 = RPAR(2) M3 = RPAR(3) M4 = RPAR(4) M5 = RPAR(5) M6 = RPAR(6) M7 = RPAR(7) XB = RPAR(8) YB = RPAR(9) XC = RPAR(10) YC = RPAR(11) C0 = RPAR(12) I1 = RPAR(13) I2 = RPAR(14) I3 = RPAR(15) I4 = RPAR(16) I5 = RPAR(17) I6 = RPAR(18) I7 = RPAR(19) D = RPAR(20) DA = RPAR(21) E = RPAR(22) EA = RPAR(23) RR = RPAR(24) RA = RPAR(25) L0 = RPAR(26) SS = RPAR(27) SA = RPAR(28) SB = RPAR(29) SC = RPAR(30) SD = RPAR(31) TA = RPAR(32) TB = RPAR(33) U = RPAR(34) UA = RPAR(35) UB = RPAR(36) ZF = RPAR(37) ZT = RPAR(38) FA = RPAR(39) MOM = RPAR(40) C SIBE = DSIN(Y(1)) SITH = DSIN(Y(2)) SIGA = DSIN(Y(3)) SIPH = DSIN(Y(4)) SIDE = DSIN(Y(5)) SIOM = DSIN(Y(6)) SIEP = DSIN(Y(7)) C COBE = DCOS(Y(1)) COTH = DCOS(Y(2)) COGA = DCOS(Y(3)) COPH = DCOS(Y(4)) CODE = DCOS(Y(5)) COOM = DCOS(Y(6)) COEP = DCOS(Y(7)) C SIBETH = DSIN(Y(1)+Y(2)) SIPHDE = DSIN(Y(4)+Y(5)) SIOMEP = DSIN(Y(6)+Y(7)) C COBETH = DCOS(Y(1)+Y(2)) COPHDE = DCOS(Y(4)+Y(5)) COOMEP = DCOS(Y(6)+Y(7)) C BEP = Y(8) THP = Y(9) PHP = Y(11) DEP = Y(12) OMP = Y(13) EPP = Y(14) C DO J = 1,7 DO I = 1,7 M(I,J) = 0D0 ENDDO ENDDO C M(1,1) = M1*RA**2 + M2*(RR**2-2*DA*RR*COTH+DA**2) + I1 + I2 M(2,1) = M2*(DA**2-DA*RR*COTH) + I2 M(2,2) = M2*DA**2 + I2 M(3,3) = M3*(SA**2+SB**2) + I3 M(4,4) = M4*(E-EA)**2 + I4 M(5,4) = M4*((E-EA)**2+ZT*(E-EA)*SIPH) + I4 M(5,5) = M4*(ZT**2+2*ZT*(E-EA)*SIPH+(E-EA)**2) + + M5*(TA**2+TB**2) + I4 + I5 M(6,6) = M6*(ZF-FA)**2 + I6 M(7,6) = M6*((ZF-FA)**2-U*(ZF-FA)*SIOM) + I6 M(7,7) = M6*((ZF-FA)**2-2*U*(ZF-FA)*SIOM+U**2) + + M7*(UA**2+UB**2) + I6 + I7 DO J = 2,7 DO I = 1,J-1 M(I,J)=M(J,I) END DO END DO C NM=7 CALL INVMAT(NM,M,MINV,EI,IP) C XD = SD*COGA + SC*SIGA + XB YD = SD*SIGA - SC*COGA + YB LANG = DSQRT ((XD-XC)**2 + (YD-YC)**2) FORCE = - C0 * (LANG - L0)/LANG FX = FORCE * (XD-XC) FY = FORCE * (YD-YC) F(1) = MOM - M2*DA*RR*THP*(THP+2*BEP)*SITH F(2) = M2*DA*RR*BEP**2*SITH F(3) = FX*(SC*COGA - SD*SIGA) + FY*(SD*COGA + SC*SIGA) F(4) = M4*ZT*(E-EA)*DEP**2*COPH F(5) = - M4*ZT*(E-EA)*PHP*(PHP+2*DEP)*COPH F(6) = - M6*U*(ZF-FA)*EPP**2*COOM F(7) = M6*U*(ZF-FA)*OMP*(OMP+2*EPP)*COOM C DO J=1,7 DO I=1,6 GP(I,J)=0D0 END DO END DO C GP(1,1) = - RR*SIBE + D*SIBETH GP(1,2) = D*SIBETH GP(1,3) = - SS*COGA GP(2,1) = RR*COBE - D*COBETH GP(2,2) = - D*COBETH GP(2,3) = - SS*SIGA GP(3,1) = - RR*SIBE + D*SIBETH GP(3,2) = D*SIBETH GP(3,4) = - E*COPHDE GP(3,5) = - E*COPHDE + ZT*SIDE GP(4,1) = RR*COBE - D*COBETH GP(4,2) = - D*COBETH GP(4,4) = - E*SIPHDE GP(4,5) = - E*SIPHDE - ZT*CODE GP(5,1) = - RR*SIBE + D*SIBETH GP(5,2) = D*SIBETH GP(5,6) = ZF*SIOMEP GP(5,7) = ZF*SIOMEP - U*COEP GP(6,1) = RR*COBE - D*COBETH GP(6,2) = - D*COBETH GP(6,6) = - ZF*COOMEP GP(6,7) = - ZF*COOMEP - U*SIEP C DO I=1,7 DY(I)=Y(I+7) ENDDO DO I=8,14 DY(I)=0.0D0 DO K=1,7 SUM=0.0D0 DO J=1,6 SUM=SUM-GP(J,K)*Y(J+14) ENDDO DY(I)=DY(I)+MINV(I-7,K)*(SUM+F(K)) ENDDO ENDDO DO I=15,20 DY(I)=0.0D0 DO J=1,7 DY(I)=DY(I)+GP(I-14,J)*Y(7+J) ENDDO ENDDO RETURN END C C SUBROUTINE JEX(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF PENDULUM PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N),RPAR(*) DO I=1,LDFY DO J=1,N DFY(I,J)=0.0D0 END DO END DO WRITE(6,*)'NO ANALYTIC COMPUTATION OF THE JACOBIAN' RETURN END C C SUBROUTINE MAS(N, AM, LMAS,RPAR,IPAR) C --- MASS MATRIX IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION AM(LMAS,N),RPAR(*) AM(1,1) = 1.0D0 AM(1,2) = 1.0D0 AM(1,3) = 1.0D0 AM(1,4) = 1.0D0 AM(1,5) = 1.0D0 AM(1,6) = 1.0D0 AM(1,7) = 1.0D0 AM(1,8) = 1.0D0 AM(1,9) = 1.0D0 AM(1,10) = 1.0D0 AM(1,11) = 1.0D0 AM(1,12) = 1.0D0 AM(1,13) = 1.0D0 AM(1,14) = 1.0D0 AM(1,15) = 0.0D0 AM(1,16) = 0.0D0 AM(1,17) = 0.0D0 AM(1,18) = 0.0D0 AM(1,19) = 0.0D0 AM(1,20) = 0.0D0 RETURN END