C----------------------------------------------------------------------- C DRIVING PROGRAM FOR THE BURGER EQUATION EXAMPLE. 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 = 2) INTEGER NPDE C NUMBER OF PDES PARAMETER (NPDE = 1) INTEGER NINTMX C MAXIMAL NUMBER OF INTEVALS ALLOWED PARAMETER (NINTMX = 2000) 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*101) 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 EXACTU(NPDE) C EXACT SOLUTION AT CERTAIN POINT 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(101) C XOUT IS A SET OF SPATIAL POINTS FOR C OUTPUT DOUBLE PRECISION COEFF C COEFF IS THE COEFFCIENT OF UXX IN THE C BURGERS' EQUATION COMMON /BURGER/ COEFF INTEGER I C----------------------------------------------------------------------- C SUBROUTINES CALLED: C BACOL C VALUES C TRUU C----------------------------------------------------------------------- C SET THE REMAINING INPUT PARAMETERS. T0 = 0.0D0 TOUT = 1.0D0 ATOL(1) = 1.D-4 RTOL(1) = ATOL(1) NINT = 10 COEFF = 1.D-3 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 MFLAG(5) = 1 C INITIALIZE IPAR and RPAR. DO I = 1, LIP IPAR(I) = 0 END DO DO I = 1, LRP RPAR(I) = 0.0D0 END DO 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, A, E8.2)') 'EPS =', COEFF, ' 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, 100 XOUT(I) = XA + DBLE(I - 1) * (XB - XA)/100.D0 30 CONTINUE XOUT(101) = XB CALL VALUES(KCOL, XOUT, NINT, X, NPDE, 101, 0, UOUT, Y, VALWRK) WRITE(6,'(/A)') 'THE OUTPUT IS ' WRITE(6,'(/A, I3, A, I4)') 'KCOL =', KCOL, ', NINT =', NINT WRITE(6,'(/A, 2A)') ' XOUT ', & ' UOUT ', ' EXACTU' DO 40 I = 1, 101 CALL TRUU(TOUT, XOUT(I), EXACTU, NPDE) WRITE(6, '(/3E18.6)') XOUT(I), UOUT(I), EXACTU(1) 40 CONTINUE GOTO 9999 100 CONTINUE WRITE(6,'(A)') 'CANNOT PROCEED DUE TO ERROR FROM BACOL.' 9999 STOP 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----------------------------------------------------------------------- DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF 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) = COEFF*UXX(1) - U(1)*UX(1) C RETURN END 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. DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF 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) = -UX(1) DFDUX(1,1) = -U(1) DFDUXX(1,1) = COEFF 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----------------------------------------------------------------------- DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF C----------------------------------------------------------------------- BVAL(1)=U(1)-0.5D0+0.5D0*TANH((-0.5D0*T-0.25D0)/(4.0D0*COEFF)) 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----------------------------------------------------------------------- DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF C----------------------------------------------------------------------- BVAL(1)=U(1)-0.5D0+0.5D0*TANH((1.D0-0.5D0*T-0.25D0)/(4.0D0*COEFF)) 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. DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF 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) = 1.0D0 DBDUX(1,1) = 0.0D0 DBDT(1) = -1.D0/COEFF/SINH(-(0.5D0*T+0.25D0)/(4.D0*COEFF))**2 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. DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF 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) = 1.0D0 DBDUX(1,1) = 0.0D0 DBDT(1) = -1.D0/COEFF/SINH(-(0.5D0*T-0.75D0)/(4.D0*COEFF))**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----------------------------------------------------------------------- DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF C----------------------------------------------------------------------- C C ASSIGN U(1:NPDE) THE INITIAL VALUES OF U(T0,X). C U(1) = 0.5D0-0.5D0*TANH((X-0.25D0)/(4.0D0*COEFF)) C RETURN END C----------------------------------------------------------------------- SUBROUTINE TRUU(T, X, U, NPDE) C----------------------------------------------------------------------- C PURPOSE: C THIS FUNCTION PROVIDES THE EXACT SOLUTION OF THE PDE. 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 C OUTPUT: DOUBLE PRECISION U(NPDE) C U(1:NPDE) IS THE EXACT SOLUTION AT THE C POINT (T,X). C----------------------------------------------------------------------- DOUBLE PRECISION COEFF COMMON /BURGER/ COEFF C----------------------------------------------------------------------- U(1) = 0.5D0-0.5D0*TANH((X-0.5D0*T-0.25D0)/(4.0D0*COEFF)) C RETURN END