C----------------------------------------------------------------------- C DRIVING PROGRAM FOR THE SYSTEM OF 4 EQUATIONS DESCRIBED IN THE PAPER. C----------------------------------------------------------------------- C CONSTANTS: INTEGER KCOL C KCOL IS THE NUMBER OF COLLOCATION POINTS C TO BE USED IN EACH SUBINTERVAL, WHICH IS C EQUAL TO THE DEGREE OF THE PIECEWISE C POLYNOMIALS MINUS ONE. C 1 < KCOL < 11. PARAMETER (KCOL = 5) INTEGER NPDE C NUMBER OF PDES PARAMETER (NPDE = 4) INTEGER NINTMX C MAXIMAL NUMBER OF INTEVALS ALLOWED PARAMETER (NINTMX = 1000) INTEGER MAXVEC C THE DIMENSION OF THE VECTOR OF C BSPLINE COEFFICIENTS PARAMETER (MAXVEC = NPDE*(NINTMX*KCOL+2)) INTEGER LRP C SEE THE COMMENT FOR RPAR PARAMETER (LRP =134+NINTMX*(35+35*KCOL+31*NPDE + +38*NPDE*KCOL+8*KCOL*KCOL)+14*KCOL + +79*NPDE+NPDE*NPDE*(21 + +4*NINTMX*KCOL*KCOL+12*NINTMX*KCOL + +6*NINTMX)) C INTEGER LIP C SEE THE COMMENT FOR IPAR PARAMETER (LIP = 115+NPDE*((2*KCOL+1)*NINTMX+4)) INTEGER LENWRK C THE DIMENSION OF ARRAY WORK WHEN WE C CALL VALUES PARAMETER (LENWRK =(KCOL+2)+KCOL*(NINTMX+1)+4) INTEGER NUOUT C THE DIMENSION OF UOUT PARAMETER (NUOUT = NPDE*21) DOUBLE PRECISION XA C THE LEFT BOUNDARY POINT PARAMETER (XA = 0.0D0) DOUBLE PRECISION XB C THE RIGHT BOUNDARY POINT PARAMETER (XB = 1.0D0) C----------------------------------------------------------------------- DOUBLE PRECISION T0 C T0 < TOUT IS THE INITIAL TIME. C DOUBLE PRECISION TOUT C TOUT IS THE DESIRED FINAL OUTPUT TIME. C DOUBLE PRECISION ATOL(NPDE) C ATOL IS THE ABSOLUTE ERROR TOLERANCE C REQUEST AND IS A SCALAR QUANTITY IF C MFLAG(2) = 0. C DOUBLE PRECISION RTOL(NPDE) C RTOL IS THE RELATIVE ERROR TOLERANCE C REQUEST AND IS A SCALAR QUANTITY IF C MFLAG(2) = 0. C INTEGER NINT C NINT IS THE NUMBER OF SUBINTERVALS C DEFINED BY THE SPATIAL MESH X. C DOUBLE PRECISION X(NINTMX+1) C X IS THE SPATIAL MESH WHICH DIVIDES THE C INTERVAL [X_A,X_B] AS: X_A = X(1) < C X(2) < X(3) < ... < X(NINT+1) = X_B. C INTEGER MFLAG(7) C THIS VECTOR OF USER INPUT DETERMINES C THE INTERACTION OF BACOL WITH DASSL. C C WORK STORAGE: DOUBLE PRECISION RPAR(LRP) C RPAR IS A FLOATING POINT WORK ARRAY C OF SIZE LRP. C INTEGER IPAR(LIP) C IPAR IS AN INTEGER WORK ARRAY C OF SIZE LIP. C C----------------------------------------------------------------------- DOUBLE PRECISION Y(MAXVEC) C ON SUCCESSFUL RETURN FROM BACOL, Y IS C THE VECTOR OF BSPLINE C COEFFICIENTS AT THE CURRENT TIME T0. C INTEGER IDID C IDID IS THE BACOL EXIT STATUS FLAG C WHICH IS BASED ON THE EXIT STATUS FROM C DASSL ON ERROR CHECKING PERFORMED BY C BACOL ON INITIALIZATION. C----------------------------------------------------------------------- DOUBLE PRECISION UOUT(NUOUT) C THE APPROXIMATION SOLUTIONS AT A SET C OF POINTS DOUBLE PRECISION VALWRK(LENWRK) C VALWRK IS A WORK ARRAY IN VALUES DOUBLE PRECISION XOUT(21) C XOUT IS A SET OF SPATIAL POINTS FOR C OUTPUT INTEGER I C----------------------------------------------------------------------- C SUBROUTINES CALLED: C BACOL C VALUES C----------------------------------------------------------------------- C SET THE REMAINING INPUT PARAMETERS. T0 = 0.0D0 TOUT = 1.8D+1 ATOL(1) = 1.D-4 RTOL(1) = ATOL(1) NINT = 10 C DEFINE THE MESH BASED ON A UNIFORM STEP SIZE. X(1) = XA DO 10 I = 2, NINT X(I) = XA + ((I-1) * (XB - XA)) / NINT 10 CONTINUE X(NINT+1) = XB C INITIALIZE THE MFLAG VECTOR. DO 20 I = 1, 7 MFLAG(I) = 0 20 CONTINUE DO I = 1, LIP IPAR(I) = 0 END DO DO I = 1, LRP RPAR(I) = 0.0D0 END DO c DO I = 2,NPDE c ATOL(I) = 1.D-4 c RTOL(I) = 1.D-4 c ENDDO WRITE(6,'(/A)') 'THE INPUT IS ' WRITE(6,'(/A, I3, A, I4, 2(A, E8.2))') 'KCOL =', KCOL, ', NINT =', & NINT, ', ATOL(1) =', ATOL(1), ', RTOL(1) =', RTOL(1) WRITE(6,'(/A, E8.2)') 'TOUT = ', TOUT CALL BACOL(T0, TOUT, ATOL, RTOL, NPDE, KCOL, NINTMX, NINT, X, & MFLAG, RPAR, LRP, IPAR, LIP, Y, IDID) C CHECK FOR AN ERROR FROM BACOL. WRITE(6,'(/A, I5)') 'IDID =', IDID IF (IDID .LT. 2) GOTO 100 XOUT(1) = XA DO 30 I = 2, 20 XOUT(I) = XA + DBLE(I - 1) * (XB - XA)/20.D0 30 CONTINUE XOUT(21) = XB CALL VALUES(KCOL, XOUT, NINT, X, NPDE, 21, 0, UOUT, Y, VALWRK) WRITE(6,'(/A)') 'THE OUTPUT IS ' WRITE(6,'(/A, I3, A, I4)') 'KCOL =', KCOL, ', NINT =', NINT WRITE(6,'(/A, 4A)') ' XOUT ', ' UOUT(1) ', & ' UOUT(2) ', ' UOUT(3) ', & ' UOUT(4)' DO 40 I = 1, 21 II = (I - 1) * 4 WRITE(6, '(/5E14.6)') XOUT(I), UOUT(II+1), UOUT(II+2), & UOUT(II+3), UOUT(II+4) 40 CONTINUE GOTO 9999 100 CONTINUE WRITE(6,'(A)') 'CANNOT PROCEED DUE TO ERROR FROM BACOL.' 9999 STOP END C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUBROUTINE DERIVF(T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THIS SUBROUTINE IS USED TO DEFINE THE INFORMATION ABOUT THE C PDE REQUIRED TO FORM THE ANALYTIC JACOBIAN MATRIX FOR THE DAE C OR ODE SYSTEM. ASSUMING THE PDE IS OF THE FORM C UT = F(T, X, U, UX, UXX) C THIS ROUTINE RETURNS THE JACOBIANS D(F)/D(U), D(F)/D(UX), AND C D(F)/D(UXX). C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION X C THE CURRENT SPATIAL COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,X). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,X). C DOUBLE PRECISION UXX(NPDE) C UXX(1:NPDE) IS THE APPROXIMATION OF THE C SECOND SPATIAL DERIVATIVE OF THE C SOLUTION AT THE POINT (T,X). C C OUTPUT: DOUBLE PRECISION DFDU(NPDE,NPDE) C DFDU(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR F C WITH RESPECT TO THE J-TH COMPONENT C OF THE UNKNOWN FUNCTION U. C DOUBLE PRECISION DFDUX(NPDE,NPDE) C DFDUX(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR F C WITH RESPECT TO THE J-TH COMPONENT C OF THE SPATIAL DERIVATIVE OF THE C UNKNOWN FUNCTION U. C DOUBLE PRECISION DFDUXX(NPDE,NPDE) C DFDUXX(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR F C WITH RESPECT TO THE J-TH COMPONENT C OF THE SECOND SPATIAL DERIVATIVE OF THE C UNKNOWN FUNCTION U. C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C ASSIGN DFDU(1:NPDE,1:NPDE), DFDUX(1:NPDE,1:NPDE), AND C DFDUXX(1:NPDE,1:NPDE) ACCORDING TO THE RIGHT HAND SIDE OF THE PDE C IN TERMS OF U(1:NPDE), UX(1:NPDE), UXX(1:NPDE). C DFDU(1,1) = - 3.D+1 * (1.D0 - U(3) - U(4)) DFDU(1,2) = 0.D0 DFDU(1,3) = 1.5D0 + 3.D+1 * U(1) DFDU(1,4) = 3.D+1 * U(1) DFDU(2,1) = 0.D0 DFDU(2,2) = - 3.D+1 * (1.D0 - U(3) - U(4)) DFDU(2,3) = 3.D+1 * U(2) DFDU(2,4) = 1.2D0 + 3.D+1 * U(2) DFDU(3,1) = 3.D+1 * (1.D0 - U(3) - U(4)) DFDU(3,2) = 0.D0 DFDU(3,3) = - 3.D+1 * U(1) - 1.5D0 - 1.D+3 * U(4) * (1.D0 - U(3) & - U(4))**2 + 2.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4)) DFDU(3,4) = - 3.D+1 * U(1) - 1.D+3 * U(3) * (1.D0 - U(3) & - U(4))**2 + 2.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4)) DFDU(4,1) = 0.D0 DFDU(4,2) = 3.D+1 * (1.D0 - U(3) - U(4)) DFDU(4,3) = - 3.D+1 * U(2) - 1.D+3 * U(3) * (1.D0 - U(3) & - U(4))**2 + 2.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4)) DFDU(4,4) = - 3.D+1 * U(2) - 1.2D0 - 1.D+3 * U(4) * (1.D0 - U(3) & - U(4))**2 + 2.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4)) C DFDUX(1,1) = - 1.D0 DFDUX(1,2) = 0.D0 DFDUX(1,3) = 0.D0 DFDUX(1,4) = 0.D0 DFDUX(2,1) = 0.D0 DFDUX(2,2) = - 1.D0 DFDUX(2,3) = 0.D0 DFDUX(2,4) = 0.D0 DFDUX(3,1) = 0.D0 DFDUX(3,2) = 0.D0 DFDUX(3,3) = 0.D0 DFDUX(3,4) = 0.D0 DFDUX(4,1) = 0.D0 DFDUX(4,2) = 0.D0 DFDUX(4,3) = 0.D0 DFDUX(4,4) = 0.D0 C DFDUXX(1,1) = 1.D-2 DFDUXX(1,2) = 0.D0 DFDUXX(1,3) = 0.D0 DFDUXX(1,4) = 0.D0 DFDUXX(2,1) = 0.D0 DFDUXX(2,2) = 1.D-2 DFDUXX(2,3) = 0.D0 DFDUXX(2,4) = 0.D0 DFDUXX(3,1) = 0.D0 DFDUXX(3,2) = 0.D0 DFDUXX(3,3) = 1.D-2 DFDUXX(3,4) = 0.D0 DFDUXX(4,1) = 0.D0 DFDUXX(4,2) = 0.D0 DFDUXX(4,3) = 0.D0 DFDUXX(4,4) = 1.D-2 C RETURN END C----------------------------------------------------------------------- SUBROUTINE DIFBXA(T, U, UX, DBDU, DBDUX, DBDT, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THE SUBROUTINE IS USED TO DEFINE THE DIFFERENTIATED BOUNDARY C CONDITIONS AT THE LEFT SPATIAL END POINT X = XA. FOR THE C BOUNDARY CONDITION EQUATION C B(T, U, UX) = 0 C THE PARTIAL DERIVATIVES DB/DU, DB/DUX, AND DB/DT ARE SUPPLIED C BY THIS ROUTINE. C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,X). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,X). C C OUTPUT: DOUBLE PRECISION DBDU(NPDE,NPDE) C DBDU(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO THE J-TH COMPONENT C OF THE UNKNOWN FUNCTION U. C DOUBLE PRECISION DBDUX(NPDE,NPDE) C DBDUX(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO THE J-TH COMPONENT C OF THE SPATIAL DERIVATIVE OF THE C UNKNOWN FUNCTION U. C DOUBLE PRECISION DBDT(NPDE) C DBDT(I) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO TIME T. C----------------------------------------------------------------------- C C ASSIGN DBDU(1:NPDE,1:NPDE), DBDU(1:NPDE,1:NPDE), AND DBDT(1:NPDE) C ACCORDING TO THE RIGHT BOUNDARY CONDITION EQUATION IN TERMS OF C U(1:NPDE), UX(1:NPDE), UXX(1:NPDE). C DBDU(1,1) = - 1.D+2 DBDU(1,2) = 0.D0 DBDU(1,3) = 0.D0 DBDU(1,4) = 0.D0 DBDU(2,1) = 0.D0 DBDU(2,2) = - 1.D+2 DBDU(2,3) = 0.D0 DBDU(2,4) = 0.D0 DBDU(3,1) = 0.D0 DBDU(3,2) = 0.D0 DBDU(3,3) = 0.D0 DBDU(3,4) = 0.D0 DBDU(4,1) = 0.D0 DBDU(4,2) = 0.D0 DBDU(4,3) = 0.D0 DBDU(4,4) = 0.D0 C DBDUX(1,1) = 1.D0 DBDUX(1,2) = 0.D0 DBDUX(1,3) = 0.D0 DBDUX(1,4) = 0.D0 DBDUX(2,1) = 0.D0 DBDUX(2,2) = 1.D0 DBDUX(2,3) = 0.D0 DBDUX(2,4) = 0.D0 DBDUX(3,1) = 0.D0 DBDUX(3,2) = 0.D0 DBDUX(3,3) = 1.D0 DBDUX(3,4) = 0.D0 DBDUX(4,1) = 0.D0 DBDUX(4,2) = 0.D0 DBDUX(4,3) = 0.D0 DBDUX(4,4) = 1.D0 C DBDT(1) = 0.D0 DBDT(2) = 0.D0 DBDT(3) = 0.D0 DBDT(4) = 0.D0 C RETURN END C----------------------------------------------------------------------- SUBROUTINE DIFBXB(T, U, UX, DBDU, DBDUX, DBDT, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THE SUBROUTINE IS USED TO DEFINE THE DIFFERENTIATED BOUNDARY C CONDITIONS AT THE RIGHT SPATIAL END POINT 1 = XB. FOR THE C BOUNDARY CONDITION EQUATION C B(T, U, UX) = 0 C THE PARTIAL DERIVATIVES DB/DU, DB/DUX, AND DB/DT ARE SUPPLIED C BY THIS ROUTINE. C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,X). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,X). C C OUTPUT: DOUBLE PRECISION DBDU(NPDE,NPDE) C DBDU(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO THE J-TH COMPONENT C OF THE UNKNOWN FUNCTION U. C DOUBLE PRECISION DBDUX(NPDE,NPDE) C DBDUX(I,J) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO THE J-TH COMPONENT C OF THE SPATIAL DERIVATIVE OF THE C UNKNOWN FUNCTION U. C DOUBLE PRECISION DBDT(NPDE) C DBDT(I) IS THE PARTIAL DERIVATIVE C OF THE I-TH COMPONENT OF THE VECTOR B C WITH RESPECT TO TIME T. C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C ASSIGN DBDU(1:NPDE,1:NPDE), DBDU(1:NPDE,1:NPDE), AND DBDT(1:NPDE) C ACCORDING TO THE RIGHT BOUNDARY CONDITION EQUATION IN TERMS OF C U(1:NPDE), UX(1:NPDE), UXX(1:NPDE). C DBDU(1,1) = 0.D0 DBDU(1,2) = 0.D0 DBDU(1,3) = 0.D0 DBDU(1,4) = 0.D0 DBDU(2,1) = 0.D0 DBDU(2,2) = 0.D0 DBDU(2,3) = 0.D0 DBDU(2,4) = 0.D0 DBDU(3,1) = 0.D0 DBDU(3,2) = 0.D0 DBDU(3,3) = 0.D0 DBDU(3,4) = 0.D0 DBDU(4,1) = 0.D0 DBDU(4,2) = 0.D0 DBDU(4,3) = 0.D0 DBDU(4,4) = 0.D0 C DBDUX(1,1) = 1.D0 DBDUX(1,2) = 0.D0 DBDUX(1,3) = 0.D0 DBDUX(1,4) = 0.D0 DBDUX(2,1) = 0.D0 DBDUX(2,2) = 1.D0 DBDUX(2,3) = 0.D0 DBDUX(2,4) = 0.D0 DBDUX(3,1) = 0.D0 DBDUX(3,2) = 0.D0 DBDUX(3,3) = 1.D0 DBDUX(3,4) = 0.D0 DBDUX(4,1) = 0.D0 DBDUX(4,2) = 0.D0 DBDUX(4,3) = 0.D0 DBDUX(4,4) = 1.D0 C DBDT(1) = 0.D0 DBDT(2) = 0.D0 DBDT(3) = 0.D0 DBDT(4) = 0.D0 C RETURN END C----------------------------------------------------------------------- SUBROUTINE BNDXA(T, U, UX, BVAL, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE C LEFT SPATIAL END POINT X = XA. C B(T, U, UX) = 0 C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,XA). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,XA). C C OUTPUT: DOUBLE PRECISION BVAL(NPDE) C BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION C AT THE LEFT BOUNDARY POINT. C----------------------------------------------------------------------- BVAL(1) = UX(1) + 1.D+2 * (2.D0-0.96D0-U(1)) BVAL(2) = UX(2) + 1.D+2 * (0.96D0-U(2)) BVAL(3) = UX(3) BVAL(4) = UX(4) C RETURN END C----------------------------------------------------------------------- SUBROUTINE BNDXB(T, U, UX, BVAL, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE C RIGHT SPATIAL END POINT X = XB. C B(T, U, UX) = 0 C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,XB). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,XB). C C OUTPUT: DOUBLE PRECISION BVAL(NPDE) C BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION C AT THE RIGHT BOUNDARY POINT. C----------------------------------------------------------------------- BVAL(1) = UX(1) BVAL(2) = UX(2) BVAL(3) = UX(3) BVAL(4) = UX(4) C RETURN END C----------------------------------------------------------------------- SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THIS SUBROUTINE DEFINES THE RIGHT HAND SIDE VECTOR OF THE C NPDE DIMENSIONAL PARABOLIC PARTIAL DIFFERENTIAL EQUATION C UT = F(T, X, U, UX, UXX). C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C DOUBLE PRECISION T C THE CURRENT TIME COORDINATE. C DOUBLE PRECISION X C THE CURRENT SPATIAL COORDINATE. C DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE APPROXIMATION OF THE C SOLUTION AT THE POINT (T,X). C DOUBLE PRECISION UX(NPDE) C UX(1:NPDE) IS THE APPROXIMATION OF THE C SPATIAL DERIVATIVE OF THE SOLUTION AT C THE POINT (T,X). C DOUBLE PRECISION UXX(NPDE) C UXX(1:NPDE) IS THE APPROXIMATION OF THE C SECOND SPATIAL DERIVATIVE OF THE C SOLUTION AT THE POINT (T,X). C C OUTPUT: DOUBLE PRECISION FVAL(NPDE) C FVAL(1:NPDE) IS THE RIGHT HAND SIDE C VECTOR F(T, X, U, UX, UXX) OF THE PDE. C----------------------------------------------------------------------- C C ASSIGN FVAL(1:NPDE) ACCORDING TO THE RIGHT HAND SIDE OF THE PDE C IN TERMS OF U(1:NPDE), UX(1:NPDE), UXX(1:NPDE). C FVAL(1) = - UX(1) + 1.D-2 * UXX(1) + 1.5D0 * U(3) - 3.D+1 * U(1) & * (1.D0 - U(3) - U(4)) FVAL(2) = - UX(2) + 1.D-2 * UXX(2) + 1.2D0 * U(4) - 3.D+1 * U(2) & * (1.D0 - U(3) - U(4)) FVAL(3) = 1.D-2 * UXX(3) + 3.D+1 * U(1) * (1.D0 - U(3) - U(4)) & - 1.5D0 * U(3) - 1.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4))**2 FVAL(4) = 1.D-2 * UXX(4) + 3.D+1 * U(2) * (1.D0 - U(3) - U(4)) & - 1.2D0 * U(4) - 1.D+3 * U(3) * U(4) * (1.D0 - U(3) & - U(4))**2 C RETURN END C----------------------------------------------------------------------- SUBROUTINE UINIT(X, U, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THIS SUBROUTINE IS USED TO RETURN THE NPDE-VECTOR OF INITIAL C CONDITIONS OF THE UNKNOWN FUNCTION AT THE INITIAL TIME T = T0 C AT THE SPATIAL COORDINATE X. C C----------------------------------------------------------------------- C SUBROUTINE PARAMETERS: C INPUT: DOUBLE PRECISION X C THE SPATIAL COORDINATE. C INTEGER NPDE C THE NUMBER OF PDES IN THE SYSTEM. C C OUTPUT: DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS VECTOR OF INITIAL VALUES OF C THE UNKNOWN FUNCTION AT T = T0 AND THE C GIVEN VALUE OF X. C----------------------------------------------------------------------- C C ASSIGN U(1:NPDE) THE INITIAL VALUES OF U(T0,X). C U(1) = 1.04D0 U(2) = 0.96D0 U(3) = 0.D0 U(4) = 0.D0 C RETURN END