#!/bin/sh # This is a shell archive (produced by GNU sharutils 4.2). # To extract the files from this archive, save it to some FILE, remove # everything before the `!/bin/sh' line above, then type `sh FILE'. # # Made on 1996-07-23 09:17 EDT by . # Source directory was `/home/wcr/final_pack'. # # Existing files will *not* be overwritten unless `-c' is specified. # This format requires very little intelligence at unshar time. # "if test", "echo", "mkdir", and "sed" may be needed. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 7362 -rw-rw-rw- p1n1.f # 7936 -rw-rw-rw- p2n1.f # 11108 -rw-rw-rw- p3n1.f # 9567 -rw-rw-rw- p1n2.f # 20697 -rw-rw-rw- p2n2.f # 9438 -rw-rw-rw- p3n2.f # 11255 -rw-rw-rw- p4n2.f # 9190 -rw-rw-rw- p1q2.f # 10718 -rw-rw-rw- p2q2.f # 8645 -rw-rw-rw- p3q2.f # 9041 -rw-rw-rw- p1q3.f # 11118 -rw-rw-rw- p2q3.f # 5997 -rw-rw-rw- p1sq1.f # 5813 -rw-rw-rw- p2sq1.f # 6997 -rw-rw-rw- p1eul3.f # 9133 -rw-rw-rw- p2eul3.f # echo=echo if mkdir _sh02162; then $echo 'x -' 'creating lock directory' else $echo 'failed to create lock directory' exit 1 fi # ============= p1n1.f ============== if test -f 'p1n1.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1n1.f' '(file already exists)' else $echo 'x -' extracting 'p1n1.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1n1.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN1 to solve the "figure 8" problem XC XC u^2 + u'^2 - 1 = 0 XC 2*u*u' - w = 0 XC XC subject to the initial conditions XC XC u(0) = 0, u'(0) = 1, w(0) = 0 XC XC and with the exact solution XC XC u(t) = sin t, w(t) = sin 2t XC XC Note that here the (KU + KW) x (KU + KW) matrix ( D_pF, D_wF ) XC has the form XC (2*p 0) XC (2*u -1) XC and hence is nonsingular only as long as p .ne. 0. This does XC not hold for t = k*pi/2. These points are removable singularities XC of the solution curve which the code passes easily but where XC DASSL fails. XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Written by W. Rheinboldt, Dec. 2 1995 XC Last modified May 25, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LX,LM,LR,LI X PARAMETER( LU=1, LW=1, LX=2*LU+LW+1, LM=LU+1 ) X PARAMETER( LR=LX*(3*LX+5)+LM*(2*LX+14), LI=2*LX ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KW = LW X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 1.0D-3 X HMIN = 0.0D0 X HMAX = 1.0D0 X NMAX = 500 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 0.0D0 X UP(1) = 1.0D0 X W(1) = 0.0D0 X T = 0.0D0 X TOUT = 8.0D0 XC XC.....Call DAEN1 XC X CALL DAEN1( FF,DFF,SOLOUT,KU,KW,U,UP,W,T,TOUT,RDATA, X & IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER ) XC XC.....Write error identifier, if needed XC X IF( IER .NE. 0 )WRITE(6,20)IER X 20 FORMAT(' Error return from DAEN1, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTN1(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,UP,W,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO X PARAMETER( ONE = 1.0D0 , TWO = 2.0D0 ) XC X FV(1) = U(1)*U(1) + UP(1)*UP(1) - ONE X FV(2) = TWO*U(1)*UP(1) - W(1) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,UP,W,T,DFV,LDF,IER ) XC X INTEGER IER,LDF X DOUBLE PRECISION U(*),UP(*),W(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO,ZER X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 , ZER = 0.0D0 ) XC X DFV(1,1) = TWO*U(1) X DFV(1,2) = TWO*UP(1) X DFV(1,3) = ZER X DFV(1,4) = ZER XC X DFV(2,1) = TWO*UP(1) X DFV(2,2) = TWO*U(1) X DFV(2,3) = -ONE X DFV(2,4) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT( TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT ) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .TRUE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE,UPE,WE,tmp XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.2D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE,UPE,WE ) XC XC.....Check for first point, write header and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),UP(1),W(1) X WRITE(6,50) UE,UPE,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),UP(1),W(1) X WRITE(6,50) UE,UPE,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF( JPOL .EQ. 0 )THEN X WRITE(6,20) NPT,T,U(1),UP(1),W(1) X WRITE(6,50) UE,UPE,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF( TASK .EQ. 'INTP' )THEN X WRITE(6,30) T,U(1),UP(1),W(1) X WRITE(6,50) UE,UPE,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),UP(1),W(1) X WRITE(6,50) UE,UPE,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEN1 for the figure-eight dae problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U= ',G18.11,' UP= ',G18.11,' W= ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U= ',G18.11,' UP= ',G18.11,' W= ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U= ',G18.11,' UP= ',G18.11,' W= ',G18.11) X 50 FORMAT(' Exact solution '/ X & ' U= ',G18.11,' UP= ',G18.11,' W= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT(T,UE,UPE,WE) XC X DOUBLE PRECISION T,UE,UPE,WE XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION SIN,COS,T2 XC XC.....Exact solution for the figure eight problem XC X T2 = T + T X UE = SIN(T) X UPE = COS(T) X WE = SIN(T2) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p1n1.f' 'failed' fi # ============= p2n1.f ============== if test -f 'p2n1.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2n1.f' '(file already exists)' else $echo 'x -' extracting 'p2n1.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2n1.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN1 to solve the (modified) XC oregenator problem XC XC a1*u1' + u1 + u1*w1 - u2 = 0 XC u2' - a2*[w1 - u2] = 0 XC u1*[1.0d0 - w1] + w1*[1.0d0 - a3*w1] = 0 XC XC with a1 = 77.27, a2 = 0.1610, a3 = 8.375D-6 subject to XC the consistent initial conditions XC XC w1(t0) = 2.70889700000D0 XC u1(t0) = w1(t0)*(1.0d0 - a3*w1(t0))/(w1(t0) - 1.0D0) XC u2(t0) = 1.000275d0 XC u1'(t0) = (1.0d0/a1)*(u2(t0) - u1(t0)*(1.0D0 + w1(t0))) XC u2'(t0) = a2*(w1(t0) - u2(t0)) XC t0 = 0.0d0 XC XC The system has a period of about T= 302.9 and is integrated from XC t = 0.0D0 to t = 305.0D0 XC XC The original Oregonator is a model of the Belousov-Zhabotinskii XC reaction proposed by Field, Koros, and Noyes, see XC XC Richard J. Field and Richard M. Noyes, "Oscillations in XC Chemical Systems IV", The J. of Chem. Phys. 60(5):1877-1884, XC March, 1974. XC XC The original model consists of a system of three ODE's which XC differs from the above system by the addition of a term eps*u3' XC on the right side of the third equation. XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by P. Isaac, April 1990 XC Last revised May 25, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LX,LM,LR,LI X PARAMETER( LU=2, LW=1, LX=2*LU+LW+1, LM=LU+1 ) X PARAMETER( LR=LX*(3*LX+5)+LM*(2*LX+14), LI=2*LX ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX,W1 XC XC.....Common array for the coefficients XC X DOUBLE PRECISION A1,A2,A3 X COMMON /COEF/ A1,A2,A3 XC XC.....Set basic dimensions XC X KU = LU X KW = LW X LRW = LR X LIW = LI XC XC.....Set coefficients XC X A1 = 77.27D0 X A2 = 0.161D0 X A3 = 8.375D-6 XC XC.....Set data arrays XC X H = 1.0D-9 X HMIN = 1.0D-10 X HMAX = 1.0D1 X NMAX = 200000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X W1 = 2.70889700000D0 X W(1) = W1 X U(1) = W1*(1.0D0 - W1*A3)/(W1 - 1.0D0) X U(2) = 1.000275D0 X UP(1) = (1.0D0/A1)*(U(2) - U(1)*(1.0D0 + W1)) X UP(2) = A2*(W1 - U(2)) X T = 0.0D0 X TOUT = 305.0D0 XC XC.....Call DAEN1 XC X CALL DAEN1( FF,DFF,SOLOUT,KU,KW,U,UP,W,T,TOUT,RDATA, X & IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER ) XC X IF( IER .NE. 0 )WRITE(6,20)IER X 20 FORMAT(' Error return from DAEN1, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTN1(6) XC X STOP X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,UP,W,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE = 1.0D0 ) XC X DOUBLE PRECISION A1,A2,A3 X COMMON /COEF/ A1,A2,A3 XC X FV(1) = A1*UP(1) + U(1)*(ONE + W(1)) - U(2) X FV(2) = UP(2) - A2*(W(1) - U(2)) X FV(3) = U(1)*(ONE - W(1)) + W(1)*(ONE - A3*W(1)) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,UP,W,T,DFV,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),UP(*),W(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO,ZER X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 , ZER = 0.0D0 ) XC X DOUBLE PRECISION A1,A2,A3 X COMMON /COEF/ A1,A2,A3 XC X DFV(1,1)= ONE + W(1) X DFV(1,2)= -ONE X DFV(1,3)= A1 X DFV(1,4)= ZER X DFV(1,5)= U(1) X DFV(1,6)= ZER XC X DFV(2,1)= ZER X DFV(2,2)= A2 X DFV(2,3)= ZER X DFV(2,4)= ONE X DFV(2,5)= -A2 X DFV(2,6)= ZER XC X DFV(3,1)= ONE - W(1) X DFV(3,2)= ZER X DFV(3,3)= ZER X DFV(3,4)= ZER X DFV(3,5)= -U(1) + ONE - TWO*A3*W(1) X DFV(3,6)= ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT( TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT ) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/1.0D0/ XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2),W(1) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF( JPOL .EQ. 0 )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF( TASK .EQ. 'INTP' )THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2),W(1) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT( ' DAEN1 for the oregonator problem' ) X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11/ X & ' UP= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11/ X & ' UP= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11/ X & ' UP= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2n1.f' 'failed' fi # ============= p3n1.f ============== if test -f 'p3n1.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p3n1.f' '(file already exists)' else $echo 'x -' extracting 'p3n1.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p3n1.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN1 to solve the batch reactor problem XC XC u1' + k2*u2*w2 = 0 XC u2' + k1*u2*u6 - km1*w4 + k2*u2*w2 = 0 XC u3' - k2*u2*w2 - k3*u4*u6 + km3*w3 = 0 XC u4' + k1*u4*u6 - km3*w3 = 0 XC u5' - k1*u2*u6 + km1*w4 = 0 XC u6' + k1*u2*u6 + k3*u4*u6 - km1*w4 - km3*w3 = 0 XC u6 - w1 + w2 + w3 + w4 - a = 0 XC w2 - (kk2*u1)/(kk2 + w1) = 0 XC w3 - (kk3*u3)/(kk3 + w1) = 0 XC w4 - (kk1*u5)/(kk1 + w1) = 0 XC XC with XC XC k1 = 21.893, km1 = 2.14D+9, k2 = 32.318, XC k3 = 21.893, km3 = 1.07D9 XC kk1 = 7.65D-18, kk2 = 4.03D-11, kk3 = 5.32D-18 XC a = 0.0131 XC XC subject to the initial conditions XC XC u1(0) = 1.5776, u2(0) = 8.32, u3(0) = 0.0, XC u4(0) = 0.0 , u5(0) = 0.0 , u6(0) = 0.0131 XC w1(0) = 7.9735D-6, w3(0) = 0.0, w4(0) = 0.0 XC XC w2(0) = (kk2*u1(0))/(kk2 + w1(0)) XC XC and the derivatives at t = 0.0 determined by the XC differential equations. XC XC The system is to be integrated from t = 0.0 to t = 1.0 [hr] XC XC The model of this batch reactor was published in XC XC M. Caracotsios and W. E. Stewart XC Sensitivity analysis of initial value problems with XC mixed ODEs and algebraic equations XC Computers and Chem. Eng. 9(4), 1985,359-365 XC XC see also XC XC L. T. Biegler and J. J. Damiano XC Nonlinear parameter estimation: A case study comparison XC AIChE Journal 32(1), 1986, 29-45 XC XC Scale factors SCM and SCT were introduced to make all variables XC dimensionless XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC W. Rheinboldt, June 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LX,LM,LR,LI X PARAMETER( LU=6, LW=4, LX=2*LU+LW+1, LM=LU+1 ) X PARAMETER( LR=LX*(3*LX+5)+LM*(2*LX+14), LI=2*LX ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER IER,JPOL,KU,KW,LIW,LRW,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Common array for the coefficients XC X DOUBLE PRECISION A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM X COMMON /COEF/A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM XC XC.......................Executable statements........................... XC XC.....Set basic dimensions XC X KU = LU X KW = LW X LRW = LR X LIW = LI XC XC.....Set the scale factors XC X SCM = 1.0D-1 X SCT = 1.0D-5 XC XC.....Set coefficients XC X A = 0.0131D0/SCM X K1 = 21.893D0*SCT*SCM X KM1 = 2.14D9*SCT X K2 = 32.318D0*SCT*SCM X K3 = 21.893D0*SCT*SCM X KM3 = 1.07D9*SCT X KK1 = 7.65D-18/SCM X KK2 = 4.03D-11/SCM X KK3 = 5.32D-18/SCM XC XC.....Set data arrays XC X H = 1.0D-9/SCT X HMIN = 1.0D-10 X HMAX = 1.0D2/SCT X NMAX = 200000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-5 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X T = 0.0D0 X TOUT = 1.0D0/SCT X U(1) = 1.5776D0/SCM X U(2) = 8.32D0/SCM X U(3) = 0.0D0 X U(4) = 0.0D0 X U(5) = 0.0D0 X U(6) = A XC X W(1) = 0.5D0*(DSQRT(( 4.0D0*U(1)+KK2)*KK2 ) - KK2) X W(2) = W(1) X W(3) = 0.0D0 X W(4) = 0.0D0 XC X UP(1) = -K2*U(2)*W(2) X UP(5) = K1*U(2)*U(6) X UP(2) = UP(1) - UP(5) X UP(3) = -UP(1) X UP(4) = 0.0D0 X UP(6) = -UP(5) XC XC.....Call DAEN1 XC X CALL DAEN1( FF,DFF,SOLOUT,KU,KW,U,UP,W,T,TOUT,RDATA, X & IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER ) XC X IF( IER .NE. 0 )WRITE(6,10)IER X 10 FORMAT(' Error return from DAEN1, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTN1(6) XC X STOP X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,UP,W,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Common array for the coefficients XC X DOUBLE PRECISION A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM X COMMON /COEF/A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM XC X FV(1) = UP(1) + K2*U(2)*W(2) X FV(2) = UP(2) + K1*U(2)*U(6) - KM1*W(4) + K2*U(2)*W(2) X FV(3) = UP(3) - K2*U(2)*W(2) - K3*U(4)*U(6) + KM3*W(3) X FV(4) = UP(4) + K1*U(4)*U(6) - KM3*W(3) X FV(5) = UP(5) - K1*U(2)*U(6) + KM1*W(4) X FV(6) = UP(6) + K1*U(2)*U(6) + K3*U(4)*U(6) - KM1*W(4) - KM3*W(3) X FV(7) = (U(6) - A) - W(1) + W(2) + W(3) + W(4) X FV(8) = W(2) - (KK2*U(1))/(KK2 + W(1)) X FV(9) = W(3) - (KK3*U(3))/(KK3 + W(1)) X FV(10) = W(4) - (KK1*U(5))/(KK1 + W(1)) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,UP,W,T,DFV,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),UP(*),W(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER I,J X DOUBLE PRECISION TMP X DOUBLE PRECISION ONE,ZER X PARAMETER( ONE = 1.0D0, ZER = 0.0D0 ) XC XC.....Common array for the coefficients XC X DOUBLE PRECISION A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM X COMMON /COEF/A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM XC X DO 20 J = 1,17 X DO 10 I = 1,10 X DFV(I,J) = ZER X 10 CONTINUE X 20 CONTINUE XC X DFV(1,2) = K2*W(2) X DFV(1,7) = ONE X DFV(1,14) = K2*U(2) XC X DFV(2,2) = K1*U(6) + K2*W(2) X DFV(2,6) = K1*U(2) X DFV(2,8) = ONE X DFV(2,14) = K2*U(2) X DFV(2,16) = -KM1 XC X DFV(3,2) = -K2*W(2) X DFV(3,4) = -K3*U(6) X DFV(3,6) = -K3*U(4) X DFV(3,9) = ONE X DFV(3,14) = -K2*U(2) X DFV(3,15) = KM3 XC X DFV(4,4) = K1*U(6) X DFV(4,6) = K1*U(4) X DFV(4,10) = ONE X DFV(4,15) = -KM3 XC X DFV(5,2) = -K1*U(6) X DFV(5,6) = -K1*U(2) X DFV(5,11) = ONE X DFV(5,16) = KM1 XC X DFV(6,2) = K1*U(6) X DFV(6,4) = K3*U(6) X DFV(6,6) = K1*U(2) + K3*U(4) X DFV(6,12) = ONE X DFV(6,15) = -KM3 X DFV(6,16) = -KM1 XC X DFV(7,6) = ONE X DFV(7,13) = -ONE X DFV(7,14) = ONE X DFV(7,15) = ONE X DFV(7,16) = ONE XC X TMP = KK2 + W(1) X DFV(8,1) = -KK2/TMP X DFV(8,13) = (KK2*U(1))/(TMP*TMP) X DFV(8,14) = ONE XC X TMP = KK3 + W(1) X DFV(9,3) = -KK3/TMP X DFV(9,13) = (KK3*U(3))/(TMP*TMP) X DFV(9,15) = ONE XC X TMP = KK1 + W(1) X DFV(10,5) = -KK1/TMP X DFV(10,13) = (KK1*U(5))/(TMP*TMP) X DFV(10,16) = ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT( TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT ) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X INTEGER I X DOUBLE PRECISION USC(6),USCP(6),WSC(4),TSC XC XC.....Common array for the coefficients XC X DOUBLE PRECISION A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM X COMMON /COEF/A,K1,KM1,K2,K3,KM3,KK1,KK2,KK3,SCT,SCM XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X TSTP = 0.2D0/SCT XC XC.....Rescale the variables XC X TSC = SCT*T X DO 10 I = 1,6 X USC(I) = U(I)*SCM X USCP(I) = UP(I)*SCM/SCT X 10 CONTINUE X DO 20 I = 1,4 X WSC(I) = W(I)*SCM X 20 CONTINUE XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,100) X WRITE(6,110) NPT,TSC X WRITE(6,130) (USC(I),I=1,6) X WRITE(6,140) (USCP(I),I=1,6) X WRITE(6,150) (WSC(I),I=1,4) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,160) TSC X WRITE(6,130) (USC(I),I=1,6) X WRITE(6,140) (USCP(I),I=1,6) X WRITE(6,150) (WSC(I),I=1,4) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF( JPOL .EQ. 0 )THEN X WRITE(6,110) NPT,TSC X WRITE(6,130) (USC(I),I=1,6) X WRITE(6,140) (USCP(I),I=1,6) X WRITE(6,150) (WSC(I),I=1,4) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF( TASK .EQ. 'INTP' )THEN X WRITE(6,120) NPT,TSC X WRITE(6,130) (USC(I),I=1,6) X WRITE(6,140) (USCP(I),I=1,6) X WRITE(6,150) (WSC(I),I=1,4) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,110) NPT,TSC X WRITE(6,130) (USC(I),I=1,6) X WRITE(6,140) (USCP(I),I=1,6) X WRITE(6,150) (WSC(I),I=1,4) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 100 FORMAT(' DAEN1 for the batch reactor problem') X 110 FORMAT(1X/' NPT = ',I6,' T= ',G18.11) X 120 FORMAT(1X/' Interpolated point',' NPT = ',I6,' T= ',G18.11) X 130 FORMAT(' U = '/(1X,4G16.8)) X 140 FORMAT(' UP = '/(1X,4G16.8)) X 150 FORMAT(' W = '/(1X,4G16.8)) X 160 FORMAT(1X/' Final point',' T= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p3n1.f' 'failed' fi # ============= p1n2.f ============== if test -f 'p1n2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1n2.f' '(file already exists)' else $echo 'x -' extracting 'p1n2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1n2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN2 to solve the index-two, XC autonomous DAE XC XC u1' = 1/u1^2 - sin(arcos u1) - 1 - w^2 XC u2' = - w XC 0 = u2 - ln u1 XC XC with the initial conditions XC XC u1(t0) = cos(0.5), u2(t0) = ln cos(0.5), w(t0) = tan(0.5), t0 = 0.5 XC XC The problem has the exact solution XC XC u1 = cos t , u2 = ln cos t, w = tan t XC XC The problem was given in XC XC R. Maerz and C. Tischendorf, Solving More General Index-2 XC Differential Equations, Computers Math. Appl. 28, 1994, 77-105 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL INIT,GF,DPGF,DWGF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=2, LW=1, LD=2 ) X PARAMETER( LR=LD*(LD+2*LU+18)+LU*(3*LU+16)+LW*(LW-9)+26 ) X PARAMETER( LI=LU+2*LD+3 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION ATOL,RTOL,RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,JNEWT,NMAX X DOUBLE PRECISION H,HMIN,HMAX XC XC.....Other parameters XC X DOUBLE PRECISION HALF X PARAMETER( HALF = 0.5D0 ) XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 0.5D-1 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 X JNEWT = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(3) = JNEWT XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set starting and end-time XC X T = HALF X TOUT = 1.5D0 XC XC.....Evaluate exact initial point XC X CALL INIT(T,U,UP,W,IER) X IF( IER .NE. 0 )THEN X WRITE(6,20)T X 20 FORMAT(' T = ',G14.6,' is an illegal starting time') X STOP X ENDIF XC X CALL DAEN2( GF,DPGF,DWGF,FF,DFF,SOLOUT, X & KU,KW,KD,U,UP,W,T,TOUT,RDATA,IDATA, X & ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,30)IER X 30 FORMAT(' Error return from DAEN2, with IER = ',I5) XC XC.....Write statistics XC X CALL WSTN2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE INIT(T,U,UP,W,IER) XC X INTEGER IER X DOUBLE PRECISION T,U(*),W(*),UP(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC INIT determines the initial point of the Maerz/Tischendorf XC index-2 problem for a given time T XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION LOG, SIN, COS, U1 XC X DOUBLE PRECISION ZER X PARAMETER( ZER = 0.0D0 ) XC X U1 = COS(T) X IF( U1 .EQ. ZER )THEN X IER = -1 X RETURN X ENDIF XC X U(1) = U1 X U(2) = LOG(U1) X UP(1) = -SIN(T) X W(1) = -UP(1)/U1 X UP(2) = -W(1) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,UP,W,T,VG,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,VG(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION SIN, ACOS, U1, ACU XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER = 0.0D0, ONE = 1.0D0 ) XC X U1 = U(1) X IF( U1 .LE. ZER )THEN X IER = -1 X RETURN X ENDIF X ACU = ACOS(U1) XC X VG(1) = (U1*U1)*UP(1) + (U1*U1)*(SIN(ACU)+W(1)*W(1)+ONE)- ONE X VG(2) = UP(2) + W(1) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DPGF(U,UP,W,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER, LDC, LDA, KA X DOUBLE PRECISION U(*),UP(*),W(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J X DOUBLE PRECISION FACT XC X FACT = U(1)*U(1) X DO 20 J = 1, KA X AC(1,J) = FACT*C(1,J) X AC(2,J) = C(2,J) X 20 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DWGF(U,UP,W,T,B,LDB,IER) XC X INTEGER IER, LDB X DOUBLE PRECISION U(*),UP(*),W(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 ) X DOUBLE PRECISION FACT XC X FACT = U(1)*U(1) XC X B(1,1) = TWO*FACT*W(1) X B(2,1) = ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER X PARAMETER( ZER = 0.0D0 ) XC X DOUBLE PRECISION LOG XC X IF( U(1) .LE. ZER )THEN X IER = -1 X RETURN X ENDIF X V(1) = U(2) - LOG(U(1)) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,A,LDA,IER) XC X INTEGER LDA, IER X DOUBLE PRECISION U(*),T,A(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER = 0.0D0, ONE = 1.0D0 ) XC X IF( U(1) .LE. ZER )THEN X IER = -1 X RETURN X ENDIF XC X A(1,1) = -ONE/U(1) X A(1,2) = ONE X A(1,3) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X INTEGER JPOL, NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT X CHARACTER*6 TASK XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION U1,U2,W1 XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Determine the exact solution XC X CALL EXCT( T,U1,U2,W1 ) XC XC.....Check for first point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20)NPT,T,U(1),U(2),W(1) X WRITE(6,50)U1,U2,W1 X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),W(1) X WRITE(6,50) U1,U2,W1 X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF( JPOL .EQ. 0 )THEN X WRITE(6,20)NPT,T,U(1),U(2),W(1) X WRITE(6,50)U1,U2,W1 X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30)T,U(1),U(2),W(1) X WRITE(6,50)U1,U2,W1 X TSEQ = TSEQ + TSTP X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20)NPT,T,U(1),U(2),W(1) X WRITE(6,50)U1,U2,W1 X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEN2 for the Maerz/Tischendorf-index-2 problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11,' W1= ',G18.11 ) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11,' W1= ',G18.11 ) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1= ',G18.11,' U2= ',G18.11,' W1= ',G18.11 ) X 50 FORMAT(' Exact solution'/ X & ' U1= ',G18.11,' U2= ',G18.11,' W1= ',G18.11 ) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,U1,U2,W1 ) XC X DOUBLE PRECISION T,U1,U2,W1 X DOUBLE PRECISION COS,LOG,TAN XC X U1 = COS(T) X U2 = LOG(U1) X W1 = TAN(T) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p1n2.f' 'failed' fi # ============= p2n2.f ============== if test -f 'p2n2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2n2.f' '(file already exists)' else $echo 'x -' extracting 'p2n2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2n2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN2 to solve the "Trajectory Prescribed XC Path Control Problem" modeled by the index 2 DAE XC XC u1' = u4*cos(u5)*sin(u6)/[R*cos(u3)] XC u2 = u4*sin(u5) XC u3' = (1/R)*u4*cos(u5)*cos(u6) XC u4' = -D(w1)/M - G*sin(u5) XC - OM^2*R*cos(u3)*(sin(u3)*cos(u5)*cos(u6) - cos(u3)*sin(u5) XC u5' = L(w1)*cos(w2)/(M*u4) + [(u4/R) - (MU/R^2 *u4)]*cos(u5) XC + 2*MU*cos(u3)*sin(u6) XC + (OM^2 *R)*[cos(u3)/u4]*[sin(u3)*cos(u6)*sin(u5) XC + cos(u3)*cos(u5)] XC u6' = L(w1)*sin(w2)/(M*u4*cos(u5)) XC + (1/R)*u4*cos(u5)*sin(u6)*tan(u3) XC - 2*MU*[cos(u3)*cos(u6)*tan(u5) - sin(u3)] XC + (OM^2 *R)*cos(u3)*sin(u3)*sin(u6)/[u4*cos(u5)] XC XC 0 = u5 + RAD*(1 + 0.0001*t^2) XC 0 = u6 - RAD*(45 + 0.001*t^2) XC XC where XC XC R = AE + u2 distance from center of earth XC AE = 20902900 radius of the earth in feet XC RHO = AR1*EXP(-AR2*u2) density of the atmosphere XC SA = 1.0 vehicle cross-sectional area XC OM = 0.72921159D-4 XC M = 2.890532728 vehicle mass XC MU = 0.1407652916D17 gravity constant XC G = MU/R^2 gravity at distance R XC RAD = pi/180 conversion factor to radians XC CL = ACL*w1/RAD aerodynamic lift coefficient XC ACL = 0.01 XC L = (1/2)*RHO*CL(w1)*SA*u4^2 XC CD = ACD1 + ACD2*CL^2 aerodynamic drag coefficient XC ACD1 = 0.04 XC ACD2 = 0.1 XC XC The initial conditions are XC XC t = 0.0 XC u1 = 0.0 (rad) , u2 = 100000 (ft), u3 = 0.0 (rad) XC u4 = 12000 (ft), u5 = -1.0 (rad) , u4 = 45.0/RAD (rad) XC XC from which the code computes XC XC w1 = 0.46650383D-01, w2= -0.91122917D-03 XC XC and XC XC u1' = 0.40394369E-03, u2' = -209.42888, u3' = 0.40394369E-03, XC u4' = -34.978260, u5' = 0.0, u6' = 0.0 XC XC The problem was given in XC XC K. E. Brenan, Stability and convergence of difference XC approximations for higher index differential algebraic XC systems with applications in trajectory control XC Ph.D. Thesis, Univ. of California at Los Angeles, 1983 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC A factor SCAL is introduced to scale the length-variables XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC.....External user subroutines XC X EXTERNAL INIT,GF,DPGF,DWGF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=6, LW=2, LD=6 ) X PARAMETER( LR=LD*(LD+2*LU+18)+LU*(3*LU+16)+LW*(LW-9)+26 ) X PARAMETER( LI=LU+2*LD+3 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION ATOL,RTOL,RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,JNEWT,NMAX X DOUBLE PRECISION H,HMIN,HMAX XC XC.....Common block for the constants XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 0.5D-1 X HMIN = 0.0D0 X HMAX = 1.0D2 X NMAX = 2000 X JPOL = 1 X JNEWT = 0 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(3) = JNEWT XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set starting and end-time XC X T = 0.0D0 X TOUT = 3.0D2 XC XC.....Set the scale factor XC X SCAL = 1.0D5 XC XC.....Evaluate initial point and set the constants XC X CALL INIT(U,UP,W,T,IER) X IF( IER .NE. 0 )THEN X WRITE(6,10) X 10 FORMAT(' Error return from the initialization') X STOP X ENDIF XC X CALL DAEN2( GF,DPGF,DWGF,FF,DFF,SOLOUT, X & KU,KW,KD,U,UP,W,T,TOUT,RDATA,IDATA, X & ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,20)IER X 20 FORMAT(' Error return from DAEN2, with IER = ',I5) XC XC.....Write statistics XC X CALL WSTN2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE INIT(U,UP,W,T,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC INIT sets the constants and the initial point for XC trajectory prescribed path control problem XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X DOUBLE PRECISION ATAN, EXP, DDIST2 XC X DOUBLE PRECISION CO5,FACT,RHO,SCAL3 X DOUBLE PRECISION H(6) XC X DOUBLE PRECISION ZER, ONE, HALF, PI, TTF1, TTF2 X PARAMETER( ZER = 0.0D0, ONE = 1.0D0, HALF = 0.5D0, X & PI = 3.141592653589793D0, X & TTF1 = 2.0D-4, TTF2 = 2.0D-3 ) XC XC.....Common block for the constants XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC XC.....Set the constants for the COEF block XC X SCAL3 = SCAL*SCAL*SCAL X ACL = 1.0D-2 X ACD1 = 4.0D-2 X ACD2 = 1.0D-1 X AE = 2.09029D7/SCAL X AMU = 0.1407653916D17/SCAL3 X AR1 = SCAL3*0.2378D-2 X AR2 = SCAL/2.38D4 X OM = 0.7292115854D-4 X RAD = PI/1.8D2 X SA = ONE/(SCAL*SCAL) X VM = 2.890532718D0 XC XC.....Set the initial U XC X U(1) = 0.0D0 X U(2) = 1.0D5/SCAL X U(3) = 0.0D0 X U(4) = 1.2D4/SCAL X U(5) = -RAD X U(6) = 4.5D1*RAD XC XC.....Set-up for the computation of W XC X UP(1) = ZER X UP(2) = ZER X UP(3) = ZER X UP(4) = ZER X UP(5) = -RAD*TTF1*T X UP(6) = RAD*TTF2*T X W(1) = ZER X W(2) = ZER XC XC.....Determine W XC X CALL GF(U,UP,W,T,H,IER) XC X W(1) = DDIST2(2,H(5),1,H(5),1,1) X W(2) = ATAN(H(6)/H(5)) XC XC.....Determine UP XC X RHO = AR1*EXP(-AR2*U(2)) X CO5 = COS(U(5)) X FACT = (VM*RAD)/(HALF*RHO*SA*U(4)*ACL) XC X UP(5) = ZER X UP(6) = ZER XC X CALL GF(U,UP,W,T,H,IER) XC X UP(1) = -H(1) X UP(2) = -H(2) X UP(3) = -H(3) X UP(4) = -H(4)/VM X UP(5) = -H(5)/FACT X UP(6) = -H(6)/(FACT*CO5) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE GF(U,UP,W,T,VG,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),VG(*),T XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC GFCT calculates the KD-dimensional vector VG = G(U,UP,W,T) XC XC Variables in the calling sequence: XC ---------------------------------- XC U D IN Array of dimension KU, the current point. XC T D IN Current time XC UP D IN Array of dimension KU, the current derivative of U XC W D IN Array of dimension KU, the current W XC VG D OUT Array of dimension KD, the vector G(U,UP,W,T) XC IER I OUT Error indicator: XC IER = 0 no error. XC IER = -1 error in GF XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION SIN, COS, EXP XC X DOUBLE PRECISION CL,CD,D,FACT,G,R,RIN,RHO,TMP,U2,U4,U4IN,W1 X DOUBLE PRECISION CO3,CO5,CO6,COW2,SI3,SI5,SI6,SIW2 XC X DOUBLE PRECISION ZER, ONE, HALF, TWO X PARAMETER( ZER = 0.0D0, ONE = 1.0D0, HALF = 0.5D0, TWO = 2.0D0 ) XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC X U2 = U(2) X U4 = U(4) X W1 = W(1) XC X R = AE + U2 X RIN = ONE/R X G = AMU*RIN*RIN X RHO = AR1*EXP(-AR2*U2) X CL = ACL*W1/RAD X CD = ACD1 + ACD2*CL*CL X D = HALF*RHO*CD*SA*U4*U4 X FACT = (VM*RAD)/(HALF*RHO*SA*U4*ACL) XC X TMP = U(3) X CO3 = COS(TMP) X SI3 = SIN(TMP) X TMP = U(5) X CO5 = COS(TMP) X SI5 = SIN(TMP) X TMP = U(6) X CO6 = COS(TMP) X SI6 = SIN(TMP) X TMP = W(2) X COW2 = COS(TMP) X SIW2 = SIN(TMP) XC X IF( (U4 .EQ. ZER) .OR. (CO3 .EQ. ZER) .OR. (CO5 .EQ. ZER) )THEN X IER = -1 X RETURN X ENDIF X U4IN = ONE/U4 XC X VG(1) = UP(1) - RIN*U4*CO5*SI6/CO3 X VG(2) = UP(2) - U4*SI5 X VG(3) = UP(3) - RIN*U4*CO5*CO6 X VG(4) = VM*UP(4) + D + VM*G*SI5 X & + OM*OM*VM*R*CO3*(SI3*CO5*CO6 - CO3*SI5) X VG(5) = FACT*UP(5) - W1*COW2 - FACT*( CO5*(U4*RIN - U4IN*G) X & + TWO*OM*CO3*SI6 X & + OM*OM*R*U4IN*CO3*(SI3*SI5*CO6 + CO3*CO5) ) X VG(6) = FACT*CO5*UP(6) - W1*SIW2 X & - FACT*( RIN*U4*CO5*CO5*SI6*SI3/CO3 X & - TWO*OM*(CO3*CO6*SI5 - SI3*CO5) X & + OM*OM*R*U4IN*CO3*SI3*SI6 ) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DPGF(U,UP,W,T,A,LDA,KA,GA,LDG,IER) XC X INTEGER IER, LDA, KA, LDG X DOUBLE PRECISION U(*),UP(*),W(*),A(LDA,*),GA(LDG,*),T XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Let DPG = DpG(U,UP,W,T) denote the KU x KU partial derivative XC of G with respect to P, and A a given KU x KA matrix XC stored in the LDA x KA array A. The routine returns the XC the product DPG*A in the LDG x MDIM array GA. XC XC Variables in the calling sequence: XC ---------------------------------- XC U D IN Array of dimension KU, the current point. XC UP D IN Array of dimension KU, the current derivative of U XC W D IN Array of dimension KU, the current W XC T D IN Current time XC A D IN Array of dimension LDA x K, the given matrix XC LDA I IN Leading dimension of A, LDA >= NVAR XC KA I IN Number of columns of A, KA >= MDIM XC GA D OUT Array of dimension KU x KA, the product DpG*A XC LDG I IN Leading dimension of GA, LDG >= MDIM XC IER I OUT Error indicator: XC ier = 0 no error. XC ier = -1 error in DPGF XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J X DOUBLE PRECISION COS,FACT,RHO XC X DOUBLE PRECISION HALF X PARAMETER( HALF = 0.5D0 ) XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC X RHO = AR1*EXP(-AR2*U(2)) X FACT = (VM*RAD)/(HALF*RHO*SA*U(4)*ACL) XC X DO 10 J = 1,KA X GA(1,J) = A(1,J) X GA(2,J) = A(2,J) X GA(3,J) = A(3,J) X GA(4,J) = VM*A(4,J) X GA(5,J) = FACT*A(5,J) X GA(6,J) = FACT*COS(U(5))*A(6,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DWGF(U,UP,W,T,A,LDA,IER) XC X INTEGER IER, LDA X DOUBLE PRECISION U(*),UP(*),W(*),A(LDA,*),T XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC DWGF evaluates the KD x KW dimensional partial derivative XC of G with respect to W and stores it in the LDA x KW XC dimensional array with LDA >= KD XC XC Variables in the calling sequence: XC ---------------------------------- XC U D IN Array of dimension KU, the current point. XC UP D IN Array of dimension KU, the current derivative of U XC W D IN Array of dimension KU, the current W XC T D IN Current time XC A D OUT Array of dimension KD x KW, the derivative DWG XC LDA D IN Leading dimension of the array A, LDA >= KD XC IER I OUT Error indicator: XC ier = 0 no error. XC ier = -1 error in DWGF. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP, COS, SIN XC X DOUBLE PRECISION CL,COW2,DCD,RHO,SIW2,TMP,U4,W1 XC X DOUBLE PRECISION ZER, HALF, TWO X PARAMETER( ZER = 0.0D0, HALF = 0.5D0, TWO = 2.0D0 ) XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC X U4 = U(4) X W1 = W(1) XC X CL = ACL*W1/RAD X DCD = TWO*ACD2*ACL*CL X RHO = AR1*EXP(-AR2*U(2)) X TMP = W(2) X COW2 = COS(TMP) X SIW2 = SIN(TMP) XC X A(1,1) = ZER X A(1,2) = ZER X A(2,1) = ZER X A(2,2) = ZER X A(3,1) = ZER X A(3,2) = ZER X A(4,1) = HALF*RHO*DCD*SA*U4*U4 X A(4,2) = ZER X A(5,1) = -COW2 X A(5,2) = W1*SIW2 X A(6,1) = -SIW2 X A(6,2) = -W1*COW2 XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,VF,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),VF(*),T XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC FF evaluates the KA = KU + KW - KD dimensional vector XC F(U,T) and stores it in the array V XC XC Variables in the calling sequence: XC ---------------------------------- XC U D IN Array of dimension KU, the current point. XC T D IN Current time XC VF D OUT Array of dimension KA containing V = F(X,T) XC IER I OUT Error indicator: XC IER = 0 no error. XC IER = -1 error in FF XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION TSQ XC X DOUBLE PRECISION ONE, FORFV, TF1, TF2 X PARAMETER( ONE = 1.0D0, FORFV = 4.5D1, X & TF1 = 1.0D-4, TF2 = 1.0D-3 ) XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC X TSQ = T*T X VF(1) = U(5)/RAD + ONE + TF1*TSQ X VF(2) = U(6)/RAD - FORFV - TF2*TSQ XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,A,LDA,IER) XC X INTEGER LDA, IER X DOUBLE PRECISION U(*),A(LDA,*),T XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC DFF evaluates the Jacobian of F at (U,T) and stores it the XC LDA x (KU+1) dimensional array A where LDA >= KA = KU+KW-KD XC XC Variables in the calling sequence: XC ---------------------------------- XC U D IN Array of dimension KU, the current point. XC T D IN Current time XC A D OUT Array of dimension KA x (KU+1) for the Jacobian XC of F at U,T. Let Fk denote the k-th component XC of F. Then the k-th row of DF should contain XC the vector of the KU+1 partial derivatives XC XC ( d/dX(1) Fk , .... , d/dX(KX) Fk, d/dT Fk ) XC XC LDA I IN Leading dimension of A, LDA >= KA XC IER I OUT Error indicator: XC ier = 0 no error. XC ier = -1 error in DFF. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER, ONE, TTF1, TTF2 X PARAMETER( ZER = 0.0D0, ONE = 1.0D0, X & TTF1 = 2.0D-4, TTF2 = 2.0D-3 ) XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC X A(1,1) = ZER X A(1,2) = ZER X A(1,3) = ZER X A(1,4) = ZER X A(1,5) = ONE/RAD X A(1,6) = ZER X A(1,7) = TTF1*T X A(2,1) = ZER X A(2,2) = ZER X A(2,3) = ZER X A(2,4) = ZER X A(2,5) = ZER X A(2,6) = ONE/RAD X A(2,7) = -TTF2*T XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X INTEGER JPOL, NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT X CHARACTER*6 TASK XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION U2,U4,UP2,UP4 XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Common block for the constants XC X DOUBLE PRECISION ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM X COMMON /COEF/ACL,ACD1,ACD2,AE,AMU,AR1,AR2, X & OM,RAD,SA,SCAL,VM XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/1.0D1/ XC XC.....Rescale the length variables XC X U2 = SCAL*U(2) X U4 = SCAL*U(4) X UP2 = SCAL*UP(2) X UP4 = SCAL*UP(4) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20)NPT,T,U(1),U2,U(3),U4,U(5),U(6),W(1),W(2), X & UP(1),UP2,UP(3),UP4,UP(5),UP(6) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U2,U(3),U4,U(5),U(6),W(1),W(2), X & UP(1),UP2,UP(3),UP4,UP(5),UP(6) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20)NPT,T,U(1),U2,U(3),U4,U(5),U(6),W(1),W(2), X & UP(1),UP2,UP(3),UP4,UP(5),UP(6) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30)T,U(1),U2,U(3),U4,U(5),U(6),W(1),W(2), X & UP(1),UP2,UP(3),UP4,UP(5),UP(6) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U2,U(3),U4,U(5),U(6),W(1),W(2), X & UP(1),UP2,UP(3),UP4,UP(5),UP(6) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEN2 for the trajectory-path-prescribed', X & ' control problem'/1X) X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G20.12,' U2 = ',G20.12,' U3 = ',G20.12/ X & ' U4 = ',G20.12,' U5 = ',G20.12,' U6 = ',G20.12/ X & ' W1 = ',G20.12,' W2 = ',G20.12/ X & ' UP1= ',G20.12,' UP2= ',G20.12,' UP3= ',G20.12/ X & ' UP4= ',G20.12,' UP5= ',G20.12,' UP6= ',G20.12) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G20.12,' U2 = ',G20.12,' U3 = ',G20.12/ X & ' U4 = ',G20.12,' U5 = ',G20.12,' U6 = ',G20.12/ X & ' W1 = ',G20.12,' W2 = ',G20.12/ X & ' UP1= ',G20.12,' UP2= ',G20.12,' UP3= ',G20.12/ X & ' UP4= ',G20.12,' UP5= ',G20.12,' UP6= ',G20.12) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G20.12,' U2 = ',G20.12,' U3 = ',G20.12/ X & ' U4 = ',G20.12,' U5 = ',G20.12,' U6 = ',G20.12/ X & ' W1 = ',G20.12,' W2 = ',G20.12/ X & ' UP1= ',G20.12,' UP2= ',G20.12,' UP3= ',G20.12/ X & ' UP4= ',G20.12,' UP5= ',G20.12,' UP6= ',G20.12) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2n2.f' 'failed' fi # ============= p3n2.f ============== if test -f 'p3n2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p3n2.f' '(file already exists)' else $echo 'x -' extracting 'p3n2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p3n2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN2 to solve the index-two problem XC XC u1' - [a - 1/(2-t)]*u1 - [(3-t)/(2-t)]*exp(t) - a*(2 - t) w = 0 XC u2' - [(a - 1)/(2-t)]*u1 + u2 - 2*exp(t) - (a - 1) w = 0 XC (2 + t)*u1 + (t^2 - 4)*u2 - (t^2 + t - 2)*exp(t) = 0 XC XC with a = 50 and the initial conditions XC XC u1(t0) = 1, u2(t0) = 1, t0 = 0, XC u1'(t0) = 1, u2'(t0) = 1, w(t0) = -1/2 XC XC The problem was originally proposed by U. Ascher and L. Petzold XC Lawrence Livermore Nat'l Lab., Num. Math. Group, Techn. XC Report UCRL-JC-107441, 1991 XC XC The exact solution is XC XC u1(t) = exp[t], u2(t) = exp[t], w = -exp[t]/(2-t) XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by W. Rheinboldt, July 1993 XC Last revised March 15, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL GF,DPGF,DWGF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=2, LW=1, LD=2 ) X PARAMETER( LR=LD*(LD+2*LU+18)+LU*(3*LU+16)+LW*(LW-9)+26 ) X PARAMETER( LI=LU+2*LD+3 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION U(LU),UP(LU),W(LW),T,TOUT XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,JNEWT,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Common array for the coefficient XC X DOUBLE PRECISION A X COMMON /COEF/ A XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set coefficient XC X A = 5.0D1 XC XC.....Set data arrays XC X H = 5.0D-3 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 X JNEWT = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(3) = JNEWT XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 1.0D0 X U(2) = 1.0D0 X UP(1) = 1.0D0 X UP(2) = 1.0D0 X W(1) = -0.5D0 X T = 0.0D0 X TOUT = 1.0D0 XC X CALL DAEN2( GF,DPGF,DWGF,FF,DFF,SOLOUT, X & KU,KW,KD,U,UP,W,T,TOUT,RDATA,IDATA, X & ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,20)IER X 20 FORMAT(' Error return from DAEN2, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTN2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,UP,W,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO,THREE X PARAMETER( ONE=1.0D0, TWO=2.0D0, THREE=3.0D0 ) XC X DOUBLE PRECISION EXP,ET,TWMT,TWMTIN XC X DOUBLE PRECISION A X COMMON /COEF/A XC X IF( T .EQ. TWO ) THEN X IER = -1 X RETURN X ENDIF XC X ET = EXP(T) X TWMT = TWO - T X TWMTIN = ONE/TWMT XC X V(1) = UP(1) - A*(U(1) + TWMT*W(1)) X & + TWMTIN*(U(1) - (THREE - T)*ET) X V(2) = UP(2) - (A - ONE)*(TWMTIN*U(1) + W(1)) + U(2) - TWO*ET XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DPGF(U,UP,W,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER, KA, LDC, LDA X DOUBLE PRECISION U(*),UP(*),W(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DWGF(U,UP,W,T,B,LDB,IER) XC X INTEGER IER, LDB X DOUBLE PRECISION U(*),UP(*),W(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO X PARAMETER( ONE=1.0D0, TWO=2.0D0 ) XC X DOUBLE PRECISION A X COMMON /COEF/A XC X B(1,1) = (T - TWO)*A X B(2,1) = ONE - A XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION TWO,FOUR X PARAMETER( TWO=2.0D0, FOUR=4.0D0 ) XC X DOUBLE PRECISION EXP,TT XC X TT = T*T X V(1) = (TWO + T)*U(1) + (TT - FOUR)*U(2) - (TT + T - TWO)*EXP(T) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,DFX,LDF,IER) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFX(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO, THREE, FOUR X PARAMETER( ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, FOUR=4.0D0 ) XC X DOUBLE PRECISION EXP,TT XC X TT = T*T X DFX(1,1) = TWO + T X DFX(1,2) = TT - FOUR X DFX(1,3) = U(1) + TWO*T*U(2) - (TT + THREE*T - ONE)*EXP(T) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X CHARACTER*6 TASK X INTEGER JPOL,NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE1,UE2,WE ) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF( TASK .EQ. 'INTP' ) THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC XC XC.....Write header XC X 10 FORMAT(' DAEN2 for the index two problem of U. Ascher') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 50 FORMAT(1X/' Exact solution'/ X & ' U1 = ',G18.11,' U2 = ',G18.11,' W= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,UE1,UE2,WE ) XC X DOUBLE PRECISION T,UE1,UE2,WE X DOUBLE PRECISION EXP, TWO X PARAMETER( TWO = 2.0D0 ) XC X UE1 = EXP(T) X UE2 = UE1 X WE = UE1/(T - TWO) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p3n2.f' 'failed' fi # ============= p4n2.f ============== if test -f 'p4n2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p4n2.f' '(file already exists)' else $echo 'x -' extracting 'p4n2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p4n2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEN2 to solve the index-two DAE XC XC u1' = -exp[t]*u1 + u2 + u4 - exp[-t] + w XC u2' = -u1 + u2 - sin[t]*u3 - cos[t] + w XC u3' = sin[t]*( u1 + u4 - sin[t] - exp[-t] ) + u3 XC u4' = cos[t]*u2 + u3 + sin[t]*u4 - exp[-t]*(1 + sin[t]) XC - cos[t]^2 - exp[t], XC 0 = u1*sin[t]^2 + u2*cos[t]^2 XC + (u3 - exp[t])*(sin[t] + 2*cos[t]) XC + sin[t]*(u4 - exp[-t])*(sin[t] + cos[t] - 1) XC - sin[t]^3 - cos[t]^3. XC XC with the initial conditions XC XC u1(t0) = 0, u2(t0) = 1, u3(t0) = 1, u4(t0) = 1, t0 = 0 XC XC u1'(t0) = 1, u2'(t0) = 0, u3'(t0) = 1, u4'(t0) = -1, w(t0) = 0 XC XC The problem has the exact solution XC XC u1 = sin[t], u2 = cos[t], u3 = exp[t], u4 = exp[-t], XC w = exp[t]*sin[t] XC XC The problem was given in XC XC K. E. Brenan XC Numerical simulation of trajectory prescribed path control XC problems by the backward differentiation method XC IEEE Trans. Auto. Control AC31, 1986, 266-269 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL GF,DPGF,DWGF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=4, LW=1, LD=4 ) X PARAMETER( LR=LD*(LD+2*LU+18)+LU*(3*LU+16)+LW*(LW-9)+26 ) X PARAMETER( LI=LU+2*LD+3 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION U(LU),UP(LU),W(LW),T,TOUT XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,JNEWT,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 0.5D-1 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 X JNEWT = 0 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(3) = JNEWT XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 0.0D0 X U(2) = 1.0D0 X U(3) = 1.0D0 X U(4) = 1.0D0 X UP(1) = 1.0D0 X UP(2) = 0.0D0 X UP(3) = 1.0D0 X UP(4) = -1.0D0 X W(1) = 0.0D0 XC X T = 0.0D0 X TOUT = 1.0D0 XC X CALL DAEN2( GF,DPGF,DWGF,FF,DFF,SOLOUT, X & KU,KW,KD,U,UP,W,T,TOUT,RDATA,IDATA, X & ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,20)IER X 20 FORMAT(' Error return from DAEN2, with IER = ',I5) XC XC.....Write statistics XC X CALL WSTN2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,UP,W,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),W(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE = 1.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET, ETIN, SIT, COT XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) XC X V(1) = UP(1) + ET*U(1) - U(2) - U(4) + ETIN - W(1) X V(2) = UP(2) + U(1) - U(2) + SIT*U(3) + COT - W(1) X V(3) = UP(3) - U(3) - SIT*( U(1) + U(4) - SIT - ETIN ) X V(4) = UP(4) - U(3) - COT*U(2) - SIT*U(4) + ETIN*( ONE + SIT ) X & + COT*COT + ET XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DPGF(U,UP,W,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER, KA, LDA, LDC X DOUBLE PRECISION U(*),UP(*),W(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J XC X IF( LDC .LT. 4 .OR. LDA .LT. 4 )THEN X IER = -1 X RETURN X ENDIF XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X AC(3,J) = C(3,J) X AC(4,J) = C(4,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DWGF(U,UP,W,T,B,LDB,IER) XC X INTEGER IER, LDB X DOUBLE PRECISION U(*),UP(*),W(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER = 0.0D0, ONE = 1.0D0 ) XC X B(1,1) = -ONE X B(2,1) = -ONE X B(3,1) = ZER X B(4,1) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET,ETIN,SIT,SIT2,COT,COT2 XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) X SIT2 = SIT*SIT X COT2 = COT*COT XC X V(1) = U(1)*SIT2 + U(2)*COT2 + ( U(3) - ET )*( SIT + TWO*COT ) X & + SIT*( U(4) - ETIN )*( SIT + COT - ONE ) X & - SIT2*SIT - COT2*COT XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,A,LDA,IER) XC X INTEGER LDA, IER X DOUBLE PRECISION U(*),T,A(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO, THREE X PARAMETER( ONE = 1.0D0, TWO = 2.0D0, THREE = 3.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET,ETIN,SIT,COT XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) XC X A(1,1) = SIT*SIT X A(1,2) = COT*COT X A(1,3) = SIT + TWO*COT X A(1,4) = SIT*( SIT + COT - ONE ) X A(1,5) = TWO*( U(1) - U(2) )*SIT*COT X & + ( U(3) - ET )*( COT - TWO*SIT ) X & - ET*( SIT + TWO*COT ) X & + COT*( U(4) - ETIN )*( SIT + COT - ONE ) X & + SIT*ETIN*( SIT + COT - ONE ) X & + SIT*( U(4) - ETIN )*( COT - SIT ) X & - THREE*SIT*COT*( SIT - COT ) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X INTEGER JPOL, NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT X CHARACTER*6 TASK XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,UE3,UE4,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE1,UE2,UE3,UE4,WE ) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20)NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEN2 for Brenans index-2 test problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11, X & ' UP3= ',G18.11,' UP4= ',G18.11/ X & ' W = ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11, X & ' UP3= ',G18.11,' UP4= ',G18.11/ X & ' W = ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11, X & ' UP3= ',G18.11,' UP4= ',G18.11/ X & ' W = ',G18.11) X 50 FORMAT(1X/' Exact solution'/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' W = ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,UE1,UE2,UE3,UE4,WE ) XC X DOUBLE PRECISION T,UE1,UE2,UE3,UE4,WE X DOUBLE PRECISION EXP, SIN, ONE X PARAMETER( ONE = 1.0D0 ) XC X UE1 = SIN(T) X UE2 = COS(T) X UE3 = EXP(T) X UE4 = ONE/UE3 X WE = UE3*UE1 XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p4n2.f' 'failed' fi # ============= p1q2.f ============== if test -f 'p1q2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1q2.f' '(file already exists)' else $echo 'x -' extracting 'p1q2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1q2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEQ2 to solve the index-two problem XC XC u1' + a*(t - 2) w = [a - 1/(2-t)]*u1 + [(3-t)/(2-t)]*exp(t) XC u2' + (1 - a ) w = [(a - 1)/(2-t)]*u1 - u2 + 2*exp(t) XC (2 + t)*u1 + (t^2 - 4)*u2 - (t^2 + t - 2)*exp(t) = 0 XC XC with a = 50 and the initial conditions XC XC u1(t0) = 1, u2(t0) = 1, t0 = 0 XC XC The problem was originally proposed by U. Ascher and L. Petzold XC Lawrence Livermore Nat'l Lab., Num. Math. Group, Techn. XC Report UCRL-JC-107441, 1991 XC XC The exact solution is XC XC u1(t) = exp[t], u2(t) = exp[t], w = exp[t]/(t-2) XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by W. Rheinboldt, July 1993 XC Last revised March 15, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,BMAT,GF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=2, LW=1, LD=2 ) X PARAMETER( LR=LU*(3*LU+15)+LD*(2*LU+17)+LW*(LW-11)+26 ) X PARAMETER( LI=2*LU+LW+2 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION U(LU),T,TOUT XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Common array for the coefficient XC X DOUBLE PRECISION A X COMMON /COEF/ A XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set coefficient XC X A = 5.0D1 XC XC.....Set data arrays XC X H = 1.0D-1 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 1.0D0 X U(2) = 1.0D0 X T = 0.0D0 X TOUT = 1.0D0 XC X CALL DAEQ2( AMAT,BMAT,GF,FF,DFF,SOLOUT,KU,KW,KD,U,T,TOUT, X & RDATA,IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,30)IER X 30 FORMAT(' Error return from DAEQ2, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTQ2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT( U,T,C,LDC,KA,AC,LDA,IER ) XC X INTEGER IER,KA,LDC,LDA X DOUBLE PRECISION U(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J XC X IF( LDC .LT. 2 .OR. LDA .LT. 2 )THEN X IER = -1 X RETURN X ENDIF XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE BMAT( U,T,B,LDB,IER ) XC X INTEGER LDB,IER X DOUBLE PRECISION U(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO X PARAMETER( ONE=1.0D0, TWO=2.0D0 ) XC X DOUBLE PRECISION A X COMMON /COEF/A XC X B(1,1) = (T - TWO)*A X B(2,1) = ONE - A XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF( U,T,V,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO,THREE X PARAMETER( ONE=1.0D0, TWO=2.0D0, THREE=3.0D0 ) XC X DOUBLE PRECISION EXP,ET,TWTIN XC X DOUBLE PRECISION A X COMMON /COEF/A XC X IF( T .EQ. TWO ) THEN X IER = -1 X RETURN X ENDIF XC X ET = EXP(T) X TWTIN = ONE/(TWO - T) XC X V(1) = (A - TWTIN)*U(1) + (THREE - T)*TWTIN*ET X V(2) = (A - ONE)*TWTIN*U(1) - U(2) + TWO*ET XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,T,V,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION TWO,FOUR X PARAMETER( TWO=2.0D0, FOUR=4.0D0 ) XC X DOUBLE PRECISION EXP,TT XC X TT = T*T X V(1) = (TWO + T)*U(1) + (TT - FOUR)*U(2) - (TT + T - TWO)*EXP(T) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,T,DFX,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFX(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO, THREE, FOUR X PARAMETER( ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, FOUR=4.0D0 ) XC X DOUBLE PRECISION EXP,TT XC X TT = T*T X DFX(1,1) = TWO + T X DFX(1,2) = TT - FOUR X DFX(1,3) = U(1) + TWO*T*U(2) - (TT + THREE*T - ONE)*EXP(T) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT( TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT ) XC X CHARACTER*6 TASK X INTEGER JPOL,NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE1,UE2,WE ) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20)NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF( TASK .EQ. 'INTP' ) THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEQ2 for the index two problem of U. Ascher') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 50 FORMAT(1X/' Exact solution'/ X & ' U1 = ',G18.11,' U2 = ',G18.11,' W= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,UE1,UE2,WE ) XC X DOUBLE PRECISION T,UE1,UE2,WE X DOUBLE PRECISION EXP, TWO X PARAMETER( TWO = 2.0D0 ) XC X UE1 = EXP(T) X UE2 = UE1 X WE = UE1/(T - TWO) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p1q2.f' 'failed' fi # ============= p2q2.f ============== if test -f 'p2q2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2q2.f' '(file already exists)' else $echo 'x -' extracting 'p2q2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2q2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEQ2 to solve the index-two, XC XC u1' - w = -exp[t]*u1 + u2 + u4 - exp[-t] XC u2' - w = -u1 + u2 - sin[t]*u3 - cos[t] XC u3' = sin[t]*( u1 + u4 - sin[t] - exp[-t] ) + u3 XC u4' = cos[t]*u2 + u3 + sin[t]*u4 - exp[-t]*(1 + sin[t]) XC - cos[t]^2 - exp[t], XC 0 = u1*sin[t]^2 + u2*cos[t]^2 XC + (u3 - exp[t])*(sin[t] + 2*cos[t]) XC + sin[t]*(u4 - exp[-t])*(sin[t] + cos[t] - 1) XC - sin[t]^3 - cos[t]^3. XC XC with the initial conditions XC XC u1(t0) = 0, u2(t0) = 1, u3(t0) = 1, u4(t0) = 1 XC XC The problem has the exact solution XC XC u1 = sin[t], u2 = cos[t], u3 = exp[t], u4 = exp[-t], XC w = exp[t]*sin[t] XC XC The problem was given in XC XC K. E. Brenan XC Numerical simulation of trajectory prescribed path control XC problems by the backward differentiation method XC IEEE Trans. Auto. Control AC31, 1986, 266-269 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,BMAT,GF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=4, LW=1, LD=4 ) X PARAMETER( LR=LU*(3*LU+15)+LD*(2*LU+17)+LW*(LW-11)+26 ) X PARAMETER( LI=2*LU+LW+2 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION U(LU),T,TOUT XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 0.5D-1 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 0.0D0 X U(2) = 1.0D0 X U(3) = 1.0D0 X U(4) = 1.0D0 X T = 0.0D0 X TOUT = 1.0D0 XC X CALL DAEQ2( AMAT,BMAT,GF,FF,DFF,SOLOUT,KU,KW,KD,U,T,TOUT, X & RDATA,IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,10)IER X 10 FORMAT(' Error return from DAEQ2, with IER = ',I5) XC XC.....Write statistics XC X CALL WSTQ2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER,KA,LDC,LDA X DOUBLE PRECISION U(*),T,C(LDC,*),AC(LDA,*) XC XC234567--x---------x---------x---------x---------x---------x---------x-- XC X INTEGER J XC X IF( LDC .LT. 4 .OR. LDA .LT. 4 )THEN X IER = -1 X RETURN X ENDIF XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X AC(3,J) = C(3,J) X AC(4,J) = C(4,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE BMAT(U,T,B,LDB,IER) XC X INTEGER LDB,IER X DOUBLE PRECISION U(*),T,B(LDB,*) XC XC234567--x---------x---------x---------x---------x---------x---------x-- XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER = 0.0D0, ONE = 1.0D0 ) XC X B(1,1) = -ONE X B(2,1) = -ONE X B(3,1) = ZER X B(4,1) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,T,VG,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,VG(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE = 1.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET, ETIN, SIT, COT XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) XC X VG(1) = -ET*U(1) + U(2) + U(4) - ETIN X VG(2) = -U(1) + U(2) - SIT*U(3) - COT X VG(3) = U(3) + SIT*( U(1) + U(4) - SIT - ETIN ) X VG(4) = U(3) + COT*U(2) + SIT*U(4) - ETIN*( ONE + SIT ) X & - COT*COT - ET XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET,ETIN,SIT,SIT2,COT,COT2 XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) X SIT2 = SIT*SIT X COT2 = COT*COT XC X V(1) = U(1)*SIT2 + U(2)*COT2 + ( U(3) - ET )*( SIT + TWO*COT ) X & + SIT*( U(4) - ETIN )*( SIT + COT - ONE ) X & - SIT2*SIT - COT2*COT XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,A,LDA,IER) XC X INTEGER LDA, IER X DOUBLE PRECISION U(*),T,A(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE, TWO, THREE X PARAMETER( ONE = 1.0D0, TWO = 2.0D0, THREE = 3.0D0 ) XC X DOUBLE PRECISION EXP, SIN, COS X DOUBLE PRECISION ET,ETIN,SIT,COT XC X ET = EXP(T) X ETIN = ONE/ET X SIT = SIN(T) X COT = COS(T) XC X A(1,1) = SIT*SIT X A(1,2) = COT*COT X A(1,3) = SIT + TWO*COT X A(1,4) = SIT*( SIT + COT - ONE ) X A(1,5) = TWO*( U(1) - U(2) )*SIT*COT X & + ( U(3) - ET )*( COT - TWO*SIT ) X & - ET*( SIT + TWO*COT ) X & + COT*( U(4) - ETIN )*( SIT + COT - ONE ) X & + SIT*ETIN*( SIT + COT - ONE ) X & + SIT*( U(4) - ETIN )*( COT - SIT ) X & - THREE*SIT*COT*( SIT - COT ) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X INTEGER JPOL, NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT X CHARACTER*6 TASK XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,UE3,UE4,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE1,UE2,UE3,UE4,WE ) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20)NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),U(3),U(4) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20)NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30)T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20)NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1) X WRITE(6,50) UE1,UE2,UE3,UE4,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEQ2 for Brenans index-2 test problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11, X & ' UP3= ',G18.11,' UP4= ',G18.11/ X & ' W = ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11, X & ' UP3= ',G18.11,' UP4= ',G18.11/ X & ' W = ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11) X 50 FORMAT(1X/' Exact solution'/ X & ' U1 = ',G18.11,' U2 = ',G18.11, X & ' U3 = ',G18.11,' U4 = ',G18.11/ X & ' W = ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,UE1,UE2,UE3,UE4,WE ) XC X DOUBLE PRECISION T,UE1,UE2,UE3,UE4,WE X DOUBLE PRECISION EXP, SIN, ONE X PARAMETER( ONE = 1.0D0 ) XC X UE1 = SIN(T) X UE2 = COS(T) X UE3 = EXP(T) X UE4 = ONE/UE3 X WE = UE3*UE1 XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--x---------x---------x---------x---------x---------x---------x-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--x---------x---------x---------x---------x---------x---------x-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2q2.f' 'failed' fi # ============= p3q2.f ============== if test -f 'p3q2.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p3q2.f' '(file already exists)' else $echo 'x -' extracting 'p3q2.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p3q2.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEQ2 to solve the index-two problem XC XC u1' + u1*w2 = u3 XC u2' + u2*w2 = u4 XC u3' + u1*w1 = 0 XC u4' + u2*w1 = -1 XC (1 - u1^2 - u2^2)/2 = 0 XC u1*u3 + u2*u4 = 0 XC XC with the initial conditions XC XC u1(0)=1, u2(0)=0, u3(0)=0, u4(0)=0, w1(0)=0, w2(0)=0 XC XC This version of a pendulum problem was originally formulated in XC XC C.W.Gear, B. Leimkuhler, G.K. Gupta XC Automatic integration of Euler-Lagrange equations with XC constraints XC J. Comp. Appl. Math. 12/13, 1985, 77-90 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by W. Rheinboldt, July 1993 XC Last revised May 25, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,BMAT,GF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=4, LW=2, LD=4 ) X PARAMETER( LR=LU*(3*LU+15)+LD*(2*LU+17)+LW*(LW-11)+26 ) X PARAMETER( LI=2*LU+LW+2 ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION U(LU),T,TOUT XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 1.0D-3 X HMIN = 0.0D0 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 1.0D0 X U(2) = 0.0D0 X U(3) = 0.0D0 X U(4) = 0.0D0 X T = 0.0D0 X TOUT = 3.55D1 XC X CALL DAEQ2( AMAT,BMAT,GF,FF,DFF,SOLOUT,KU,KW,KD,U,T,TOUT, X & RDATA,IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,30)IER X 30 FORMAT(' Error return from DAEQ2, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTQ2(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER,KA,LDC,LDA X DOUBLE PRECISION U(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X AC(3,J) = C(3,J) X AC(4,J) = C(4,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE BMAT(U,T,B,LDB,IER) XC X INTEGER LDB,IER X DOUBLE PRECISION U(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER X PARAMETER( ZER=0.0D0 ) XC X B(1,1) = ZER X B(2,1) = ZER X B(3,1) = U(1) X B(4,1) = U(2) XC X B(1,2) = U(1) X B(2,2) = U(2) X B(3,2) = ZER X B(4,2) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER=0.0D0, ONE=1.0D0 ) XC X V(1) = U(3) X V(2) = U(4) X V(3) = ZER X V(4) = -ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION HALF X PARAMETER( HALF=0.5D0 ) XC X V(1) = HALF - HALF*(U(1)*U(1) + U(2)*U(2)) X V(2) = U(1)*U(3) + U(2)*U(4) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,DFX,LDF,IER) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFX(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER X PARAMETER( ZER=0.0D0 ) XC X DFX(1,1) = -U(1) X DFX(1,2) = -U(2) X DFX(1,3) = ZER X DFX(1,4) = ZER X DFX(1,5) = ZER X DFX(2,1) = U(3) X DFX(2,2) = U(4) X DFX(2,3) = U(1) X DFX(2,4) = U(2) X DFX(2,5) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X CHARACTER*6 TASK X INTEGER JPOL,NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.5D0/ XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1),W(2) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1),W(2) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF( JPOL .EQ. 0 )THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1),W(2) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1),W(2) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4),W(1),W(2) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEQ2 for an index two pendulum problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G16.9/ X & ' U1= ',G16.9,' U2= ',G16.9, X & ' U3= ',G16.9,' U4= ',G16.9/ X & ' UP1= ',G16.9,' UP2= ',G16.9, X & ' UP3= ',G16.9,' UP4= ',G16.9/ X & ' W1= ',G16.9,' W2= ',G16.9) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1= ',G16.9,' U2= ',G16.9, X & ' U3= ',G16.9,' U4= ',G16.9/ X & ' UP1= ',G16.9,' UP2= ',G16.9, X & ' UP3= ',G16.9,' UP4= ',G16.9/ X & ' W1= ',G16.9,' W2= ',G16.9) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1= ',G16.9,' U2= ',G16.9, X & ' U3= ',G16.9,' U4= ',G16.9/ X & ' UP1= ',G16.9,' UP2= ',G16.9, X & ' UP3= ',G16.9,' UP4= ',G16.9/ X & ' W1= ',G16.9,' W2= ',G16.9) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p3q2.f' 'failed' fi # ============= p1q3.f ============== if test -f 'p1q3.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1q3.f' '(file already exists)' else $echo 'x -' extracting 'p1q3.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1q3.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEQ3 to solve the linear test problem XC XC ( 2 -1 ) (u1") + ( 1 ) w = ( 0 ) XC ( 0 1 ) (u2") + ( -1) ( -t) XC t u1 + u2 = 1 XC XC subject to the initial conditions XC XC u1(0) = 1, u2(0) = 1, u1'(0) = 0, u2'(0) = -1. XC XC which has the exact solution XC XC u1(t) = -(1/12)t^3 + 1, u1'(t) = -(1/4)t^2 XC u2(t) = (1/12)t^4 - t + 1, u2'(t) = (1/3)t^3 - 1 XC w(t) = t^2 + t XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,BMAT,GF,FF,DFF,D2FF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=2, LW=1, LD=2 ) X PARAMETER( LR = LU*(4*LU+19)+LD*(LD+4*LU+34)-LW*(2*LU+28)+43 ) X PARAMETER( LI=2*(LU+2) ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 1.0D-1 X HMIN = 1.0D-12 X HMAX = 0.5D0 X NMAX = 10000 X JPOL = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 1.0D0 X U(2) = 1.0D0 X UP(1) = 0.0D0 X UP(2) = -1.0D0 X T = 0.0D0 X TOUT = 6.0D0 XC XC.....Call DAEUL3 XC X CALL DAEQ3( AMAT,BMAT,GF,FF,DFF,D2FF,SOLOUT,KU,KW,KD, X & U,UP,W,T,TOUT,RDATA,IDATA,ATOL,RTOL, X & RWORK,LRW,IWORK,LIW,IER) XC XC.....Write error identifier, if needed XC X IF(IER .NE. 0) WRITE(6,30)IER X 30 FORMAT(' Error return from DAEQ3, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTQ3(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,UP,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER,KA,LDC,LDA X DOUBLE PRECISION U(*),UP(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J X DOUBLE PRECISION TWO X PARAMETER( TWO=2.0D0 ) XC X DO 10 J = 1, KA X AC(1,J) = TWO*C(1,J) - C(2,J) X AC(2,J) = C(2,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE BMAT(U,UP,T,B,LDB,IER) XC X INTEGER LDB,IER X DOUBLE PRECISION U(*),UP(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE=1.0D0 ) XC X B(1,1) = ONE X B(2,1) = -ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,UP,T,VG,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),T,VG(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER X PARAMETER( ZER=0.0D0 ) XC X VG(1) = ZER X VG(2) = -T XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE=1.0D0 ) XC X FV(1) = T*U(1) + U(2) - ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,T,DFV,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE=1.0D0 ) XC X DFV(1,1) = T X DFV(1,2) = ONE X DFV(1,3) = U(1) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE D2FF( U,T,V,D2FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),D2FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION TWO X PARAMETER( TWO=2.0D0 ) XC X D2FV(1) = TWO*V(1)*V(3) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC NOTE: This sample routine assumes increasing times. For XC decreasing times the step TSTP has to be negative and the XC test for interpolation has to be modified accordingly XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .TRUE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,UPE1,UPE2,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.2D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT(T,UE1,UE2,UPE1,UPE2,WE) XC XC.....Print the start point and initialize TSEQ XC X IF (TASK .EQ. 'START') THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEQ3 for a linear test problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G16.8/ X & ' X = ',2G18.10/' DX = ',2G18.10,' W = ',G18.9) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.10/ X & ' X = ',2G18.10/' DX = ',2G18.10' W = ',G18.9) X 40 FORMAT(1X/' Final point',' T= ',G18.10/ X & ' X = ',2G18.10/' DX = ',2G18.10' W = ',G18.9) X 50 FORMAT(' Exact solution'/ X & ' X = ',2G18.10/' DX = ',2G18.10' W = ',G18.9) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT(T,UE1,UE2,UPE1,UPE2,WE) XC X DOUBLE PRECISION T,T2,T3,UE1,UE2,UPE1,UPE2,WE XC X DOUBLE PRECISION ONE,A1,A2,A3 X PARAMETER( ONE=1.0D0 ) X PARAMETER( A1=1.0D0/1.2D1, A2=1.0D0/4.0D0, A3=1.0D0/3.0D0 ) XC X T2 = T*T X T3 = T2*T X UE1 = -A1*T3 + ONE X UE2 = A1*T2*T2 - T + ONE X UPE1= -A2*T2 X UPE2= A3*T3 - ONE X WE = T2 + T XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p1q3.f' 'failed' fi # ============= p2q3.f ============== if test -f 'p2q3.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2q3.f' '(file already exists)' else $echo 'x -' extracting 'p2q3.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2q3.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEQ3 to solve the nonlinear test problem XC XC u1" + u1*w = exp[t]*(1 + sin[t]) XC u2" + u2*w = 1/(1+t)*(2/(1+t)^2 + sin[t]) XC 0 = exp[t]/(1+t)^3 - 3*exp[t]/(1+t)^2 + exp[1/(1+t)]/(1+t) XC -(u1*u2^2 - 3*u1*u2 + exp[u2])*u2 XC XC with the initial conditions at t0 = 0 XC XC u1(t0) = 1, u2(t0) = 1, u1'(t0) = 1, u2'(t0) = -1, XC XC The problem has the exact solution XC XC u1 = exp[t], u2 = (1+t)^(-1), w = sin[t] XC XC The DAE is a modified version of a problem given in XC XC U.M. Ascher, and L.R. Petzold, Projected Implicit Runge-Kutta XC Methods for Differential Algebraic Equations, XC Lawrence Livermore National Lab., Numerical Mathematics Group, XC Technical Report UCRL-JC-104037 XC XC The modification was constructed so that the existence condition XC XC ( A(U,P,T) B(U,P,T) ) XC rank ( ) = KU + KW XC ( D_uF(U,P,T) 0 ) XC XC is violated at u1 = exp[s], u2 = (1+s)^(-1) with XC s .= 0.0375627.. Thus, when started at t = 0 the process XC fails near t = s. When started beyond the singular point XC there are no difficulties. XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,BMAT,GF,FF,DFF,D2FF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LW,LD,LR,LI X PARAMETER( LU=2, LW=1, LD=2 ) X PARAMETER( LR = LU*(4*LU+19)+LD*(LD+4*LU+34)-LW*(2*LU+28)+43 ) X PARAMETER( LI=2*(LU+2) ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION ATOL,RTOL,RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU),W(LW) XC XC.....Other variables XC X INTEGER KD,KU,KW,LIW,LRW,IER,JPOL,NMAX X DOUBLE PRECISION H,HMIN,HMAX X DOUBLE PRECISION EXP XC XC.....Set basic dimensions XC X KU = LU X KW = LW X KD = LD X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 0.1D-1 X HMIN = 1.0D-12 X HMAX = 1.0D1 X NMAX = 2000 X JPOL = 0 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X T = 0.05D0 X TOUT = 2.0D0 X U(1) = EXP(T) X U(2) = 1.0D0/(1.0D0+T) X UP(1) = U(1) X UP(2) = -U(2)*U(2) XC X CALL DAEQ3( AMAT,BMAT,GF,FF,DFF,D2FF,SOLOUT,KU,KW,KD, X & U,UP,W,T,TOUT,RDATA,IDATA,ATOL,RTOL, X & RWORK,LRW,IWORK,LIW,IER) XC X IF(IER .NE. 0) WRITE(6,30)IER X 30 FORMAT(' Error return from DAEQ3, with IER = ',I5) XC XC.....Write statistics XC X CALL WSTQ3(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,UP,T,C,LDC,KA,AC,LDA,IER) XC X INTEGER IER,KA,LDC,LDA X DOUBLE PRECISION U(*),UP(*),T,C(LDC,*),AC(LDA,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X INTEGER J XC X DO 10 J = 1, KA X AC(1,J) = C(1,J) X AC(2,J) = C(2,J) X 10 CONTINUE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE BMAT(U,UP,T,B,LDB,IER) XC X INTEGER LDB,IER X DOUBLE PRECISION U(*),UP(*),T,B(LDB,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X B(1,1) = U(1) X B(2,1) = U(2) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,UP,T,VG,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),T,VG(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP,SIN,SIT,OT,OTI XC X DOUBLE PRECISION ZER,ONE,TWO X PARAMETER( ZER = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) XC X OT = ONE + T X IF( OT .NE. ZER )THEN X OTI = ONE/OT X ELSE X IER = -1 X RETURN X ENDIF X SIT = SIN(T) XC X VG(1) = EXP(T)*(ONE + SIT) X VG(2) = (TWO*OTI*OTI + SIT)*OTI XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,T,FV,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP X DOUBLE PRECISION OT,OTI,U2 XC X DOUBLE PRECISION ZER, ONE, THREE X PARAMETER( ZER = 0.0D0, ONE = 1.0D0, THREE = 3.0D0 ) XC X OT = ONE + T X IF( OT .NE. ZER )THEN X OTI = ONE/OT X ELSE X IER = -1 X RETURN X ENDIF X U2 = U(2) XC X FV(1) = OTI*( EXP(T)*OTI*(OTI - THREE) + EXP(OTI) ) X & - U2*( U(1)*U2*(U2 - THREE) + EXP(U2) ) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,T,DFV,LDF,IER) XC X INTEGER LDF, IER X DOUBLE PRECISION U(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP X DOUBLE PRECISION OT,OTI,U2 XC X DOUBLE PRECISION ZER,ONE,TWO,THREE,SEVEN X PARAMETER( ZER=0.0D0, ONE = 1.0D0, TWO=2.0D0 ) X PARAMETER( THREE = 3.0D0, SEVEN =7.0D0 ) XC X OT = ONE + T X IF( OT .NE. ZER )THEN X OTI = ONE/OT X ELSE X IER = -1 X RETURN X ENDIF X U2 = U(2) XC X DFV(1,1) = U2*U2*(THREE - U2) X DFV(1,2) = THREE*U(1)*U2*(TWO - U2) - (ONE + U2)*EXP(U2) X DFV(1,3) = -OTI*OTI*( EXP(T)*(OTI*(THREE*OTI - SEVEN) + THREE) X & + EXP(OTI)*(OTI + ONE) ) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE D2FF(U,T,V,D2FV,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),D2FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP X DOUBLE PRECISION OT,OTI,U2 XC X DOUBLE PRECISION ZER,ONE,TWO,THREE,FOUR,SIX,TWLVE,THRTEN X PARAMETER( ZER=0.0D0, ONE = 1.0D0, TWO=2.0D0, THREE=3.0D0 ) X PARAMETER( FOUR=4.0D0, SIX = 6.0D0, TWLVE =1.2D1, THRTEN=1.3D1 ) XC X OT = ONE + T X IF( OT .NE. ZER )THEN X OTI = ONE/OT X ELSE X IER = -1 X RETURN X ENDIF X U2 = U(2) XC X D2FV(1) = V(1)*V(2)*SIX*U2*(TWO - U2) X & + V(2)*V(2)*( SIX*U(1)*(ONE-U2) - (TWO+U2)*EXP(U2) ) X & + V(3)*V(3)*OTI*OTI* X & ( EXP(T)*(OTI*(TWLVE*OTI*(OTI-TWO)+THRTEN)-THREE) X & + EXP(OTI)*OTI*(OTI*(OTI+FOUR)+TWO) ) X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,W,T,TLAST,TNEXT) XC X INTEGER JPOL, NPT X DOUBLE PRECISION U(*),UP(*),W(*),T,TLAST,TNEXT X CHARACTER*6 TASK XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Local variables XC X DOUBLE PRECISION UE1,UE2,UPE1,UPE2,WE XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Evaluate the exact solution XC X CALL EXCT( T,UE1,UE2,UPE1,UPE2,WE ) XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20)NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50)UE1,UE2,UPE1,UPE2,WE X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2),W(1) X WRITE(6,50) UE1,UE2,UPE1,UPE2,WE X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEQ3 for a modified Ascher-Petzold problem') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 40 FORMAT(1X/' Final point',' T= ',G18.11/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) X 50 FORMAT(1X/' Exact solution'/ X & ' U1 = ',G18.11,' U2 = ',G18.11/ X & ' UP1= ',G18.11,' UP2= ',G18.11,' W= ',G18.11) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT( T,UE1,UE2,UPE1,UPE2,WE ) XC X DOUBLE PRECISION T,UE1,UE2,UPE1,UPE2,WE XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP, SIN, ONE X PARAMETER( ONE = 1.0D0 ) XC X UE1 = EXP(T) X UE2 = ONE/(ONE + T) X UPE1 = UE1 X UPE2 = -UE2*UE2 X WE = SIN(T) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7==C XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2q3.f' 'failed' fi # ============= p1sq1.f ============== if test -f 'p1sq1.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1sq1.f' '(file already exists)' else $echo 'x -' extracting 'p1sq1.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1sq1.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAESQ1 to solve the simple test problem XC XC u1' = u2 XC u3' = 1 XC 0 = u1 - u3*u3 + 2*u2 XC XC subject to the initial conditions XC XC u1(0) = 11, u2(0) = -5, u3(0) = 1 XC XC The exact solution is XC XC u1(t) = 6 exp[(1-t)/2] + t^2 - 4*t + 8 XC u2(t) = -3 exp[(1-t)/2] + 2*t - 4 XC u3(t) = t XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by W. Rheinboldt, April 1992 XC Last revised March 14, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,GF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LA,LR,LI X PARAMETER( LU=3, LA=1 ) X PARAMETER( LR=LU*(3*LU+15)-9*LA, LI=3*LU ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION U(LU),RWORK(LR),RDATA(10) XC XC.....Other variables XC X INTEGER IER,IMPAS,ITARG,LIW,LRW,KU,KA,NMAX X DOUBLE PRECISION UTARG X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KA = LA X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 1.0D-3 X HMIN = 0.0D0 X HMAX = 1.0D0 X UTARG = 3.0D0 X NMAX = 2000 X ITARG = 3 X IMPAS = 0 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X RDATA(4) = UTARG X IDATA(1) = NMAX X IDATA(3) = ITARG X IDATA(4) = IMPAS XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U(1) = 11.0D0 X U(2) = -5.0D0 X U(3) = 1.0D0 XC XC.....Call DAESQ1 XC X CALL DAESQ1( AMAT,GF,FF,DFF,SOLOUT,KU,KA,U,RDATA, X & IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER ) XC X IF( IER .NE. 0 )WRITE(6,30)IER X 30 FORMAT(' Error return from DAESQ1, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTSQ1(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,UP,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X V(1) = UP(1) X V(2) = UP(3) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE X PARAMETER( ONE = 1.0D0 ) X X V(1) = U(2) X V(2) = ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION TWO X PARAMETER( TWO = 2.0D0 ) XC X V(1) = U(1) - U(3)*U(3) + TWO*U(2) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,DFU,LDF,IER) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),DFU(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ONE,TWO X PARAMETER( ONE = 1.0D0, TWO = 2.0D0 ) XC X DFU(1,1) = ONE X DFU(1,2) = TWO X DFU(1,3) = -TWO*U(3) XC X IER = 0 X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,NPT,U,GAMMA) XC X CHARACTER*6 TASK X INTEGER NPT X DOUBLE PRECISION U(*),GAMMA XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION U3,UE1,UE2 XC XC.....Evaluate the exact solution XC X U3 = U(3) X CALL EXCT( U3,UE1,UE2 ) XC XC.....Check for first point, write header XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,U(1),U(2),U(3) X WRITE(6,40) UE1,UE2,U3 X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,30) U(1),U(2),U(3) X WRITE(6,40) UE1,UE2,U3 X RETURN X ENDIF XC XC.....Write out the computed point XC X WRITE(6,20) NPT,U(1),U(2),U(3) X WRITE(6,40) UE1,UE2,U3 XC XC.....Formats XC X 10 FORMAT(' DAESQ1 for a simple test problem'/1X) X 20 FORMAT(' NPT = ',I6,' U ='/(1X,3G18.11)) X 30 FORMAT(1X/' Final point, U= '/(1X,3G18.11)) X 40 FORMAT(' Exact solution '/1X,3G18.11) XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE EXCT(U3,UE1,UE2) XC X DOUBLE PRECISION U3,UE1,UE2 XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION EXP,TMP XC XC.....Exact solution for the simple autonomous problems XC X TMP = 0.5D0*(1.0D0-U3) X UE1 = 6.0D0*EXP(TMP) + ((U3-4.0D0)*U3 + 8.0D0) X UE2 = -3.0D0*EXP(TMP) + 2.0D0*U3 - 4.0D0 XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p1sq1.f' 'failed' fi # ============= p2sq1.f ============== if test -f 'p2sq1.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2sq1.f' '(file already exists)' else $echo 'x -' extracting 'p2sq1.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2sq1.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAESQ1 to solve the autonomous DAE XC XC u3' = u4 XC u4' = u2 XC 0 = u1 + u2 + u3 XC 0 = u4 - u1*u1 - cc XC XC with variable parameter cc and initial conditions XC XC The problem was originally considered in XC XC F. Takens, Constrained equations: A study of implicit XC differential equations and their discontinuous solutions, XC Lecture Notes in Math. Vol 525, Springer Verlag 1976, pp143-234 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Originally written by W. Rheinboldt, April 1992 XC Last revised March 14, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL AMAT,GF,FF,DFF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LA,LD,LR,LI X PARAMETER( LU=4, LA=2, LD=2 ) X PARAMETER( LR=LU*(3*LU+15)-9*LA+3*LD+(LD+2)*(LD+2), LI=3*LU ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION U(LU),RWORK(LR),RDATA(10) XC XC.....Other variables XC X INTEGER IER,IMPAS,ITARG,LIW,LRW,KU,KA,NMAX X DOUBLE PRECISION U1,U2,UTARG X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Common array for the coefficient XC X DOUBLE PRECISION CC X COMMON /COEF/CC XC XC.....Set basic dimensions XC X KU = LU X KA = LA X LRW = LR X LIW = LI XC XC.....Set coefficient XC X CC = 1.0D0 XC XC.....Set tolerance and step data XC X H = 1.0D-2 X HMIN = 0.0D0 X HMAX = 1.0D0 X UTARG = 3.0D0 X NMAX = 250 X ITARG = 3 X IMPAS = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X RDATA(4) = UTARG X IDATA(1) = NMAX X IDATA(3) = ITARG X IDATA(4) = IMPAS XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X U1 = -0.2D0 X U2 = 1.0D0 X U(1) = U1 X U(2) = U2 X U(3) = -U1 - U2 X U(4) = CC + U1*U1 XC XC.....Call DAESQ1 XC X CALL DAESQ1( AMAT,GF,FF,DFF,SOLOUT,KU,KA,U,RDATA, X & IDATA,ATOL,RTOL,RWORK,LRW,IWORK,LIW,IER ) XC X IF(IER .NE. 0) WRITE(6,10)IER X 10 FORMAT(' Error return from DAESQ1, with ier = ',I5) XC X CALL WSTSQ1(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE AMAT(U,UP,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X V(1) = UP(3) X V(2) = UP(4) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF(U,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X V(1) = -U(4) X V(2) = -U(2) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF(U,V,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),V(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION CC X COMMON /COEF/CC XC X V(1) = U(1) + U(2) + U(3) X V(2) = U(4) - U(1)*U(1) - CC XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF(U,DFU,LDF,IER) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),DFU(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER,ONE,TWO X PARAMETER( ZER=0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) XC X DFU(1,1) = ONE X DFU(1,2) = ONE X DFU(1,3) = ONE X DFU(1,4) = ZER X DFU(2,1) = -TWO*U(1) X DFU(2,2) = ZER X DFU(2,3) = ZER X DFU(2,4) = ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,NPT,U,GAMMA) XC X CHARACTER*6 TASK X INTEGER NPT X DOUBLE PRECISION U(*),GAMMA XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION CC X COMMON /COEF/CC XC XC.....Check for first point, write header XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10)CC X WRITE(6,20) NPT,U(1),U(2),U(3),U(4) X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) U(1),U(2),U(3),U(4) X RETURN X ENDIF XC XC.....Write out the computed point XC X WRITE(6,20) NPT,U(1),U(2),U(3),U(4) X WRITE(6,30)GAMMA XC XC.....Formats XC X 10 FORMAT(' DAESQ1 for the Takens problem '/ X & ' with CC = ',G14.6/1X) X 20 FORMAT(' NPT = ',I6,' U ='/(1X,4G18.10)) X 30 FORMAT(' GAMMA= ',G18.11) X 40 FORMAT(1X/' Final point, U= '/(1X,4G18.10)) XC X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2sq1.f' 'failed' fi # ============= p1eul3.f ============== if test -f 'p1eul3.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p1eul3.f' '(file already exists)' else $echo 'x -' extracting 'p1eul3.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p1eul3.f' && XC234567==1=========2=========3=========4=========5=========6=========7= XC XC Sample driver for using DAEUL3 to solve the autonomous XC pendulum problem XC XC u1" + 2u1 w = 0 XC u2" + 2u2 w = -1 XC u1^2 + u2^2 - 1 = 0 XC XC subject to the initial conditions XC XC u1(t0) = 0, u2(t0) = 1, t0 = 0, u1'(t0) = 1, u2'(t0) = 0. XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC.....External user subroutines XC X EXTERNAL MMAT,GF,FF,DFF,D2FF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LA,LR,LI X PARAMETER( LU=2, LA=1 ) X PARAMETER( LR=LU*(4*LU+48) - 28*LA + 44 ) X PARAMETER( LI=2*(LU+2) ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU) XC XC.....Other variables XC X INTEGER KU,KA,LIW,LRW,IER,JPOL,MID,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Set basic dimensions XC X KU = LU X KA = LA X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 5.0D-2 X HMIN = 1.0D-12 X HMAX = 0.5D0 X NMAX = 10000 X JPOL = 1 X MID = 1 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(4) = MID XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set initial point XC X T = 0.0D0 X U(1) = 0.0D0 X U(2) = -1.0D0 X UP(1) = 1.0D0 X UP(2) = 0.0D0 X TOUT = 6.0D0 XC XC.....Call DAEUL3 XC X CALL DAEUL3( MMAT,GF,FF,DFF,D2FF,SOLOUT,KU,KA, X & U,UP,T,TOUT,RDATA,IDATA,ATOL,RTOL, X & RWORK,LRW,IWORK,LIW,IER ) XC XC.....Write error identifier, if needed XC X IF(IER .NE. 0) WRITE(6,10)IER X 10 FORMAT(' Error return from DAEUL3, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTEUL(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE MMAT(U,T,V,MV,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),MV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X MV(1) = V(1) X MV(2) = V(2) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE GF( U,UP,T,GV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),T,GV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X DOUBLE PRECISION ZER, ONE X PARAMETER( ZER=0.0D0, ONE=1.0D0 ) XC X GV(1) = ZER X GV(2) = -ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE FF( U,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X DOUBLE PRECISION ONE X PARAMETER( ONE=1.0D0 ) XC X FV(1) = U(1)*U(1) + U(2)*U(2) - ONE XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE DFF( U,T,DFV,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X DOUBLE PRECISION TWO,ZER X PARAMETER( TWO=2.0D0, ZER=0.0D0 ) XC X DFV(1,1) = TWO*U(1) X DFV(1,2) = TWO*U(2) X DFV(1,3) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE D2FF( U,T,V,D2FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),D2FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X DOUBLE PRECISION TWO X PARAMETER( TWO=2.0D0 ) XC X D2FV(1) = TWO*(V(1)*V(1) + V(2)*V(2)) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,T,TLAST,TNEXT) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .TRUE. ) XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points. XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.1D0/ XC XC.....Print the start point and initialize TSEQ XC X IF( TASK .EQ. 'START' )THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),UP(1),UP(2) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),UP(1),UP(2) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),UP(1),UP(2) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEUL3 for the pendulum') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G18.10/ X & ' U = ',2G18.10/' UP = ',2G18.10) X 30 FORMAT(1X/' Interpolated point',' T= ',G18.10/ X & ' U = ',2G18.10/' UP = ',2G18.10) X 40 FORMAT(1X/' Final point',' T= ',G18.10/ X & ' U = ',2G18.10/' UP = ',2G18.10) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7= XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7= SHAR_EOF : || $echo 'restore of' 'p1eul3.f' 'failed' fi # ============= p2eul3.f ============== if test -f 'p2eul3.f' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'p2eul3.f' '(file already exists)' else $echo 'x -' extracting 'p2eul3.f' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'p2eul3.f' && XC234567==1=========2=========3=========4=========5=========6=========7== XC XC Sample driver for using DAEUL3 to solve the following nonautonomous XC double-pendulum problem with rotating pivot-point XC (in dimensionless form) XC XC u1'' + 2*[(u1 - b sin t) - (u3 - u1)]*w = 0 XC u2'' + 2*[(u2 - b cos t) - (u4 - u2)]*w = g XC u3'' + 2*qm*(u3 - u1)*w = 0 XC u4'' + 2*qm*(u4 - u2)*w = g XC (u1 - sin t)^2 + (u2 - cos t)^2 = r1^2 XC (u3 - u1)^2 + (u4 - u2)^2 = r2^2 XC XC subject to the initial conditions XC XC u1(0) = 0, u2(0) = 1 + r1, u3(0) = 0, u4(0) = 1 + r1 + r2 XC u1'(0) = 0, u2'(0) = 0, u3'(0) = 0, u4'(0) = 0, XC XC Here g = 5, qm = 0.5, b = 0.1, r1 = 4, r2 = 8 XC XC Link with DAEPAK, MANPAK and MANAUX XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC Written by W. Rheinboldt, Feb. 23, 1996 XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....External user subroutines XC X EXTERNAL MMAT,GF,FF,DFF,D2FF,SOLOUT XC XC.....Set array dimensions XC X INTEGER LU,LA,LR,LI X PARAMETER( LU=4, LA=2 ) X PARAMETER( LR=LU*(4*LU+48) - 28*LA + 44 ) X PARAMETER( LI=2*(LU+2) ) X INTEGER IWORK(LI),IDATA(10) X DOUBLE PRECISION RWORK(LR),RDATA(10) X DOUBLE PRECISION T,TOUT,U(LU),UP(LU) XC XC.....Other variables XC X INTEGER KU,KA,LIW,LRW,IER,JPOL,MID,NMAX X DOUBLE PRECISION ATOL,RTOL,H,HMIN,HMAX XC XC.....Common array for the constants XC X DOUBLE PRECISION G,B,R1,R2 X COMMON /COEF/ G,B,R1,R2 XC XC.....Set basic dimensions XC X KU = LU X KA = LA X LRW = LR X LIW = LI XC XC.....Set data arrays XC X H = 1.0D-3 X HMIN = 1.0D-12 X HMAX = 1.0D0 X NMAX = 10000 X JPOL = 1 X MID = 0 XC X RDATA(1) = H X RDATA(2) = HMIN X RDATA(3) = HMAX X IDATA(1) = NMAX X IDATA(2) = JPOL X IDATA(4) = MID XC XC.....Set tolerances XC X RTOL = 1.0D-8 X ATOL = RTOL*1.0D-1 XC XC.....Set constants XC X G = 5.0D0 X B = 0.1D0 X R1 = 4.0D0 X R2 = 8.0D0 XC XC.....Set initial point XC X T = 0.0D0 X U(1) = 0.0D0 X U(2) = 1.0D0+R1 X U(3) = 0.0D0 X U(4) = 1.0D0+R1+R2 X UP(1) = 0.0D0 X UP(2) = 0.0D0 X UP(3) = 0.0D0 X UP(4) = 0.0D0 X TOUT = 1.2D1 XC XC.....Call DAEUL3 XC X CALL DAEUL3( MMAT,GF,FF,DFF,D2FF,SOLOUT,KU,KA, X & U,UP,T,TOUT,RDATA,IDATA,ATOL,RTOL, X & RWORK,LRW,IWORK,LIW,IER ) XC XC.....Write error identifier, if needed XC X IF(IER .NE. 0) WRITE(6,50)IER X 50 FORMAT(' Error return from DAEUL3, with ier = ',I5) XC XC.....Write statistics XC X CALL WSTEUL(6) XC X STOP XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE MMAT(U,T,V,MV,IER) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),MV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X MV(1) = V(1) X MV(2) = V(2) X MV(3) = V(3) X MV(4) = V(4) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE GF( U,UP,T,GV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),UP(*),T,GV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION ZER X PARAMETER( ZER=0.0D0 ) XC X DOUBLE PRECISION G,B,R1,R2 X COMMON /COEF/ G,B,R1,R2 XC X GV(1) = ZER X GV(2) = G X GV(3) = ZER X GV(4) = G XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE FF( U,T,FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION COS,SIN XC X DOUBLE PRECISION G,B,R1,R2 X COMMON /COEF/ G,B,R1,R2 XC X FV(1) = (U(1) - SIN(T))**2 + (U(2) - COS(T))**2 - R1*R1 X FV(2) = (U(3) - U(1))**2 + (U(4) - U(2))**2 - R2*R2 XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE DFF( U,T,DFV,LDF,IER ) XC X INTEGER LDF,IER X DOUBLE PRECISION U(*),T,DFV(LDF,*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION COS,SIN,CT,ST XC X DOUBLE PRECISION ZER,TWO X PARAMETER( ZER=0.0D0, TWO=2.0D0 ) XC X DOUBLE PRECISION G,B,R1,R2 X COMMON /COEF/ G,B,R1,R2 XC X CT = B*COS(T) X ST = B*SIN(T) X DFV(1,1) = TWO*(U(1) - ST) X DFV(1,2) = TWO*(U(2) - CT) X DFV(1,3) = ZER X DFV(1,4) = ZER X DFV(1,5) = -TWO*(U(1)*CT - U(2)*ST) XC X DFV(2,1) = -TWO*(U(3) - U(1)) X DFV(2,2) = -TWO*(U(4) - U(2)) X DFV(2,3) = TWO*(U(3) - U(1)) X DFV(2,4) = TWO*(U(4) - U(2)) X DFV(2,5) = ZER XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE D2FF( U,T,V,D2FV,IER ) XC X INTEGER IER X DOUBLE PRECISION U(*),T,V(*),D2FV(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X DOUBLE PRECISION COS,SIN,CT,ST,TMP XC X DOUBLE PRECISION TWO,FOUR X PARAMETER( TWO=2.0D0, FOUR=4.0D0 ) XC X DOUBLE PRECISION G,B,R1,R2 X COMMON /COEF/ G,B,R1,R2 XC X CT = B*COS(T) X ST = B*SIN(T) X TMP = TWO*(V(1)*V(1) + V(2)*V(2)) X D2FV(1) = TMP - FOUR*(CT*V(1)*V(5) - ST*V(2)*V(5)) X & + TWO*(U(1)*ST + U(2)*CT)*V(5)*V(5) X D2FV(2) = TMP + TWO*(V(3)*V(3) + V(4)*V(4)) X & - FOUR*(V(1)*V(3) + V(2)*V(4)) XC X IER = 0 X RETURN XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE SOLOUT(TASK,JPOL,NPT,U,UP,T,TLAST,TNEXT) XC X CHARACTER*6 TASK X INTEGER NPT,JPOL X DOUBLE PRECISION U(*),UP(*),T,TLAST,TNEXT XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC.....Parameter for intermediate prints between interpolated points XC X LOGICAL BETWN X PARAMETER( BETWN = .FALSE. ) XC XC.....Control variables to be saved XC X DOUBLE PRECISION TSEQ,TSTP X SAVE TSEQ,TSTP XC XC.....Step between interpolated points XC.....Note that for decreasing times TSTP should be negative XC X DATA TSTP/0.2D0/ XC XC.....Print the start point and initialize TSEQ XC X IF (TASK .EQ. 'START') THEN X WRITE(6,10) X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4) X TSEQ = T + TSTP X RETURN X ENDIF XC XC.....Check for final point XC X IF( TASK .EQ. 'FINAL' )THEN X WRITE(6,40) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4) X RETURN X ENDIF XC XC.....Check for JPOL = 0 when each point is to be printed XC X IF (JPOL .EQ. 0) THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4) X RETURN X ENDIF XC XC.....Check if an interpolated point is to be printed or if XC.....a point has been computed with T between TLAST and TSEQ XC X IF (TASK .EQ. 'INTP') THEN X WRITE(6,30) T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4) X ELSEIF( TASK .EQ. 'PRNT' .AND. BETWN )THEN X IF( (TLAST .LE. T .AND. T .LE. TSEQ) .OR. X & (TLAST .GE. T .AND. T .GE. TSEQ) )THEN X WRITE(6,20) NPT,T,U(1),U(2),U(3),U(4), X & UP(1),UP(2),UP(3),UP(4) X ENDIF X ENDIF XC XC.....Test for interpolation XC X IF( (TLAST .LE. TSEQ .AND. TSEQ .LE. TNEXT) .OR. X & (TLAST .GE. TSEQ .AND. TSEQ .GE. TNEXT) )THEN X T = TSEQ X TASK = 'INTP' X TSEQ = TSEQ + TSTP X ELSE X TASK = 'PRNT' X ENDIF XC X RETURN XC XC.....Formats XC X 10 FORMAT(' DAEUL3 for a wobbling double pendulum') X 20 FORMAT(1X/' NPT = ',I6,' T= ',G16.8/ X & ' X = ',4G16.8/ X & ' DX = ',4G16.8) X 30 FORMAT(1X/' Interpolated point',' T= ',G16.8/ X & ' X = ',4G16.8/ X & ' DX = ',4G16.8) X 40 FORMAT(1X/' Final point',' T= ',G16.8/ X & ' X = ',4G16.8/ X & ' DX = ',4G16.8) XC X END XC XC234567==1=========2=========3=========4=========5=========6=========7== XC X SUBROUTINE ERROUT(KL, LOUT) XC X INTEGER KL, LOUT(*) XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC XC Variables in the calling sequence XC --------------------------------- XC KL I OUT Number of different output units to be used XC by MSGPRT. For KL <= 0 and KL > 5 all printout XC by MSGPRT is suppressed. XC LOUT I OUT Array of dimension KL, 1 <= KL<= 5, for the XC KL output-units to be used by MSGPRT. XC XC234567--1---------2---------3---------4---------5---------6---------7-- XC X KL = 1 X LOUT(1) = 6 XC XC.....End of LUNIT XC X RETURN X END XC XC234567==1=========2=========3=========4=========5=========6=========7== SHAR_EOF : || $echo 'restore of' 'p2eul3.f' 'failed' fi rm -fr _sh02162 exit 0