SUBROUTINE INITAL C... C... THE DAVIDENKO DIFFERENTIAL EQUATION CAN BE WRITTEN IN MATRIX C... FORM AS C... C... - - - - - C... J*DX/DT = -F, X(T0) = X0 (1)(2) C... C... WHERE C... C... - C... J JACOBIAN MATRIX FOR THE NONLINEAR EQUATIONS C... (NLE) C... - C... X SOLUTION VECTOR FOR THE DAVIDENKO ODE C... C... T INDEPENDENT VARIABLE IN THE DAVIDENKO METHOD C... C... T0 INITIAL VALUE OF T C... C... - C... X0 INITIAL CONDITION VECTOR FOR EQUATION (1) C... C... - - C... IN THE CODING IN SUBROUTINE DERV, J IS STORED IN ARRAY A, AND F IS C... STORED IN ARRAY B. THE UNCOUPLING OF THE DERIVATIVES IN EQUATION C... (1) IS DONE IN SUBROUTINE DERV BY A CALL TO SUBROUTINES DECOMP AND C... SOLVE (LISTED AND DISCUSSED IN DETAIL IN FORSYTHE, G. E., M. A. C... MALCOLM AND C. B. MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPU- C... TATAIONS, CHAPTER 3, PRENTICE-HALL, ENGLEWOOD-CLIFFS, N.J., 1977). C... ANY OTHER ROUTINE(S) FOR LINEAR ALGEBRAIC EQUATIONS CAN BE USED IN C... -- - C... PLACE OF DECOMP AND SOLVE (THE GENERAL PROBLEM IS AX = B). C... C... - - C... THE PROGRAMMING OF EQUATION (1) IN TERMS OF VECTORS X, DX/DT AND C... - - C... F, AND MATRIX J, IS INTENDED TO ILLUSTRATE THE GENERAL FEATURES OF C... DAVIDENKO*S METHOD, I.E., THESE VECTORS AND MATRIX CAN BE EXPANDED C... TO ACCOMMODATE A PROBLEM OF ANY SIZE. HERE WE DIMENSION FOR C... TWO NLE PROVIDED BY L. T. BIEGLER: C... C... X1**2 + X2**2 = 17 C... C... 2*X1**(1/3) + X2**(1/2) = 4 C... C... DOUBLE PRECISION IS USED TO POSSIBLY ENHANCE THE PERFORMANCE C... OF DAVIDENKO*S METHOD WHEN THE JACOBIAN MATRIX IS POORLY C... CONDITIONED IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(N=2) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(N) 2 /F/ DXDT(N) 3 /CV/ A(N,N), B(N), WORK(N), IPVT(N), 4 COND C... C... INITIAL ESTIMATE OF THE SOLUTION IF(NORUN.EQ.1)THEN X(1)=1.0D0 X(2)=1.0D0 ELSE IF(NORUN.EQ.2)THEN X(1)=0.1D0 X(2)=5.0D0 END IF C... C... INITIALIZE ALL OF THE MODEL CALCULATIONS CALL DERV RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(N=2) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(N) 2 /F/ DXDT(N) 3 /CV/ A(N,N), B(N), WORK(N), IPVT(N), 4 COND C... C... JACOBIAN MATRIX A(1,1)=2.0D0*X(1) A(1,2)=2.0D0*X(2) A(2,1)=2.0D0*(1.0D0/3.0D0)*X(1)**(-2.0D0/3.0D0) A(2,2)= (1.0D0/2.0D0)*X(2)**(-1.0D0/2.0D0) C... C... COMPUTE THE RIGHT HAND SIDE VECTOR FOR THE DAVIDENKO DIFFERENTIAL C... EQUATION. NOTE THAT THESE STATEMENTS CONSTITUTE A DEFINITION OF C... THE ELEMENTS OF B (OR VECTOR F IN EQUATION (1)). THEY FOLLOW C... DIRECTLY FROM THE STATEMENT OF THE NONLINEAR EQUATIONS B(1)=-( X(1)**2 +X(2)**2 -17.0D0) B(2)=-(2.0D0*X(1)**(1.0D0/3.0D0)+X(2)**(1.0D0/2.0D0)- 4.0D0) C... C... SOLVE FOR THE VECTOR OF DERIVATIVES, DX/DT, IN EQUATION (1) (IN C... EFFECT, EQUATION (1) IS PREMULTIPLIED BY THE INVERSE OF JACOBIAN C... MATRIX J) CALL DECOMP(N,N,A,COND,IPVT,WORK) C... C... CHECK THE CONDITION OF MATRIX A BEFORE COMPUTING THE DERIVATIVE C... VECTOR, DX/DT. IF A IS NEAR SINGULAR, TERMINATE THE CURRENT RUN IF(COND.LT.1.0D+20)GO TO 3 WRITE(NO,4)COND 4 FORMAT(1H ,//, 1 30H THE CONDITION NUMBER OF A IS ,E10.3,33H SO THE CURRENT RUN IS 1 TERMINATED) NSTOP=1 C... C... COMPUTE THE DERIVATIVE VECTOR DX/DT (INDIRECTLY VIA ARRAY B) 3 CALL SOLVE(N,N,A,B,IPVT) C... C... TRANSFER ARRAY B TO ARRAY DXDT SO THE DERIVATIVE VECTOR, DX/DT, C... IS AVAILABLE IN COMMON/F/ DO 2 I=1,N DXDT(I)=B(I) 2 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(N=2) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(N) 2 /F/ DXDT(N) 3 /CV/ A(N,N), B(N), WORK(N), IPVT(N), 4 COND C... C... DIMENSION THE ARRAYS FOR PLOTTING THE EQUATION RESIDUALS AND C... CONDITION NUMBER DIMENSION F1R(11), F2R(11), CONP(11), TP(11) C... C... MONITOR PROGRESS OF THE CALCULATION WRITE(*,*)' T = ',T C... C... INITIALIZE OUTPUT IF(T.LT.0.001D0)THEN C... C... PRINT HEADING FOR THE SOLUTION WRITE(NO,1) 1 FORMAT(9X, 'T',8X,'X1', 8X,'X2',/, + 10X, 4X,'DX1/DT',4X,'DX2/DT',/, + 10X, 6X,'RES1', 6X,'RES2',/, + 10X, 5X,'CONDITION NUMBER',/) C... C... INITIALIZE COUNTER FOR PLOTTED SOLUTION IP=0 END IF C... C... PRINT THE NUMERICAL SOLUTION. THE ELEMENTS OF B MUST FIRST BE C... RECOMPUTED SINCE B CONTAINS THE VECTOR DX/DT AFTER THE CALL TO C... SOLVE IN SUBROUTINE DERV 3 B(1)=-( X(1)**2 +X(2)**2 -17.0D0) B(2)=-(2.0D0*X(1)**(1.0D0/3.0D0)+X(2)**(1.0D0/2.0D0)- 4.0D0) C... C... PRINT THE SOLUTION, DERIVATIVES OF EQUATION (1) AND THE RESI- C... DUALS OF THE NONLINEAR EQUATIONS AT THE BEGINNING AND END OF C... EACH RUN IP=IP+1 IF((IP.EQ.1).OR.(IP.EQ.11))THEN WRITE(NO,2)T,(X(I),I=1,N),(DXDT(I),I=1,N),(B(I),I=1,N),COND 2 FORMAT(F10.1,2F10.5,/,10X,2F10.5,/,10X,2F10.5,/,16X,E11.3,/) END IF C... C... STORE THE RESIDUALS AND CONDITION NUMBER FOR PLOTTING F1R(IP)=DLOG10(DABS(B(1))) F2R(IP)=DLOG10(DABS(B(2))) CONP(IP)=DLOG10(DABS(COND)) TP(IP)=T C... C... PLOT THE LOG OF THE RESIDUALS OF THE NLE TO INVESTIGATE THE C... EXPONENTIAL CONVERGENCE OF DAVIDENKO*S ALGORITHM IF(IP.LT.11)RETURN CALL SPLOTS(1,IP,TP,F1R) WRITE(NO,11) CALL SPLOTS(1,IP,TP,F2R) WRITE(NO,12) CALL SPLOTS(1,IP,TP,CONP) WRITE(NO,13) 11 FORMAT(1H ,//,16X,'LOG10 OF RESIDUAL OF EQUATION 1 VS T') 12 FORMAT(1H ,//,16X,'LOG10 OF RESIDUAL OF EQUATION 2 VS T') 13 FORMAT(1H ,//,16X,'LOG10 OF CONDITION NUMBER VS T') RETURN END