c c revision date: August 30, 1993 c c----------------------------------------------------------------------- c *** author *** c c author and contact c c luca dieci (*) c school of mathematics c georgia institute of technology c atlanta, ga 30332 c phone: (404) 853-9209 c e-mail: dieci@math.gatech.edu c c (*) please communicate promptly c any mistake you should find. c c----------------------------------------------------------------------- c c this file contains two sample drivers, driver1 and driver2 c to be used for solving differential riccati equations c with the solver DRESOL (available from netlib). c (users must compile one driver at the time) c c -- driver 1 begins -- program driver1 c This is a sample driver for the routine DRESOL. c It sets up the basic calling sequence for DRESOL and c solves a stiff (2,2) differential Riccati equation c (it is Example 4 in the paper "Numerical Integration c of the Differential Riccati Equation and some Related Issues", c SIAM J. Numer. Analysis 29-3 (1992), pp. 781-815). implicit real (a-h,o-z) external PROBL dimension NEQ(6),X(2,2),Y(36),RWORK(114),IWORK(24) c This common block is in DRESOL. We need it, because we want c to know the value of TCJ, which is where the most current Jacobian c factorization occurred. common /INTDRE/ REUVCT, REUVPT, TCJ, TPJ, NNQ, NNP, ICONTR, 1 NONLIN, NSUPER, NN, NNS, NQS, NPS, NEQN1, NEQN11, NEQN12, 2 NEQN21, NEQN22, NP3, NP23, NP23Q, NTT, NTTP, NTT2P, NTT1, 3 NTT1Q, NTT2Q, NFLEVC, NFLEVP c dimension of A11 NNQ=2 c dimension of A22 NNP=2 NEQN=NNP*NNQ c it is a stiff problem, use the BDF with Newton's like iteration MF=21 c it is an unsymmetric DRE ICONTR=1 c want to perform quasi-Newton iteration NONLIN=1 c want superstability check because problem has turning point potential NSUPER=1 c NEQ(1)=NEQN NEQ(2)=NNQ NEQ(3)=NNP NEQ(4)=ICONTR NEQ(5)=NONLIN NEQ(6)=NSUPER c give initial time, T, final time, TOUT, and initial conditions T=-1.e0 TOUT=1.e0 do 3 j=1,2 do 3 i=1,2 3 X(i,j)=0.e0 c set up dimensions of arrays NX=2 NY=36 NEQLEN=6 LREX=114 LIW=24 c give mixed error tolerances ATOL=1.e-4 RTOL=ATOL c first call to the solver ISTATE=1 c go up to TOUT ITASK=1 c no additional input given IOPT=0 c ready to call the solver CALL DRESOL(NEQ,NEQLEN,X,NX,Y,NY,MF,T,TOUT,RTOL,ATOL,ITASK, * ISTATE,IOPT,RWORK,LREX,IWORK,LIW,PROBL,RARR,IARR) c if ISTATE is negative, we stop the iteration if(ISTATE.le.0)then write(6,*)' ISTATE is ', ISTATE stop endif c for comparison purposes, for t>0, the exact solution is c c | t/2 | sqrt(epsilon) | c X(1)= | -|- | (however, there are other -unstable- c | 0 | sqrt(epsilon) | solution trajectories). c c print the solution at T=TOUT write(6,43) 43 format(2x,'Current T and current solution are '/' ') write(6,50)t write(6,50)((X(i,j),j=1,nnq),i=1,nnp) 50 format(2x,2e22.10) c Since we have used MF=21, we have the eigenvalues of the blocks A11 and c A22 of the transformed system (not the most recent ones typically, c but those computed at the predictor for the quasi-Newton iteration). c Check if we found purely imaginary or unstable eigenvalues write(6,51)IWORK(1),IWORK(2) 51 format(/ /,2x,'IWORK(1) : ',i4,2x,' IWORK(2) : ',i4) c Print eigenvalue information. N1=3+NNP+2*NNP*NNP+2*NNQ*NNQ N2=N1+NNP N3=N2+NNP+NNQ N4=N3+NNQ c the value of LWM is LWM = (MAXORD+1)*NEQN+20 if(MF.eq.21.or.MF.eq.20)MAXRD=6 if(MF.eq.11.or.MF.eq.10)MAXRD=13 c check against having set a non-default max-order if(IOPT.eq.1.and.IWORK(5).ne.0)MAXRD=IWORK(5)+1 LWM=MAXRD*NEQN+20 c print the value of t where eigenvalues have been last computed write(6,52)TCJ 52 format(/ /,2x,'Last eigenvalues computation at TCJ: ',e22.10) write(6,53) 53 format(/ /,2x,'Eigenvalues of A11(t)+A12(t)*X(t) '/' ') write(6,50)(RWORK(LWM+N3+I-1),RWORK(LWM+N4+I-1),I=1,NNQ) write(6,54) 54 format(/ /,2x,'Eigenvalues of A22(t)-X(t)*A12(t) '/' ') write(6,50)(RWORK(LWM+N1+I-1),RWORK(LWM+N2+I-1),I=1,NNP) c print final statistics on the run write(6,55) 55 format(/ /,2x,'Final statistics on this run '/' ') write(6,60)IWORK(11),IWORK(12),IWORK(13),RWORK(11),RWORK(12) 60 format 1 (2x,'nst= ',i5,' nfe=',i5,' nje=',i5,' hu=',e9.3,' hcur=',e9.3) stop end c subroutine PROBL(T,A,RARR,IARR) implicit real(a-h,o-z) dimension A(4,4) epsi=1.e-5 NN=4 do 10 j=1,NN do 10 i=1,NN A(i,j)=0.0e0 10 continue onep=1.e0/epsi A(1,1)=-t*onep*.5e0 A(3,1)=0.5e0 A(3,2)=1.e0 A(4,2)=1.e0 A(1,3)=onep A(2,4)=onep A(3,4)=t*onep*.5e0 return end c c running this program produces the following output c on a Sun Sparc station with machine precision about 1.e-7 c c Current T and current solution are c c 0.1000000000E+01 c 0.5000000000E+00 0.3162277397E-02 c -0.5743791402E-14 0.3162277630E-02 c c c IWORK(1) : 0 IWORK(2) : 0 c c c Last eigenvalues computation at TCJ: 0.2070293665E+01 c c c Eigenvalues of A11(t)+A12(t)*X(t) c c 0.0000000000E+00 0.0000000000E+00 c 0.3162278748E+03 0.0000000000E+00 c c c Eigenvalues of A22(t)-X(t)*A12(t) c c -0.1035146797E+06 0.0000000000E+00 c -0.3162241211E+03 0.0000000000E+00 c c c Final statistics on this run c c nst= 74 nfe= 141 nje= 85 hu=0.165E+01 hcur=0.165E+01 c c -- driver 1 ended -- c c -- driver 2 begins -- program driver2 c This is a sample driver for the routine DRESOL. c It sets up the basic calling sequence for DRESOL and c solves a symmetric (6,6) DRE (an optimal control problem) c (it is Example 5 in the paper "Numerical Integration c of the Differential Riccati Equation and some Related Issues", c SIAM J. Numer. Analysis 29-3 (1992), pp. 781-815). implicit real (a-h,o-z) external PROBL dimension NEQ(6),X(6,6),Y(324),RWORK(526),IWORK(32) common /INTDRE/ REUVCT, REUVPT, TCJ, TPJ, NNQ, NNP, ICONTR, 1 NONLIN, NSUPER, NN, NNS, NQS, NPS, NEQN1, NEQN11, NEQN12, 2 NEQN21, NEQN22, NP3, NP23, NP23Q, NTT, NTTP, NTT2P, NTT1, 3 NTT1Q, NTT2Q, NFLEVC, NFLEVP c dimension of A11 NNQ=6 c dimension of A22 NNP=6 NEQN=NNP*NNQ c use the BDF with Newton's like iteration MF=21 c it is a symmetric DRE ICONTR=0 c want to perform the chord iteration NONLIN=0 c no need for superstability check because constant coefficients NSUPER=0 c NEQ(1)=NEQN NEQ(2)=NNQ NEQ(3)=NNP NEQ(4)=ICONTR NEQ(5)=NONLIN NEQ(6)=NSUPER c give initial time, T, final time, TOUT, and initial conditions T=0.e0 TOUT=1.e10 do 3 j=1,6 do 3 i=1,6 3 X(i,j)=0.e0 c set up dimensions of arrays NX=6 NY=324 NEQLEN=6 LREX=526 LIW=32 c give mixed error tolerances ATOL=1.e-1 RTOL=ATOL c first call to the solver ISTATE=1 c go up to TOUT ITASK=1 c no additional input given IOPT=0 c ready to call the solver CALL DRESOL(NEQ,NEQLEN,X,NX,Y,NY,MF,T,TOUT,RTOL,ATOL,ITASK, * ISTATE,IOPT,RWORK,LREX,IWORK,LIW,PROBL,RARR,IARR) c if ISTATE is negative, we stop the iteration if(ISTATE.le.0)then write(6,*)' ISTATE is ', ISTATE stop endif c print the solution at T=TOUT write(6,43) 43 format(2x,'Current T and current solution are '/' ') write(6,50)t write(6,50)((X(i,j),j=1,nnq),i=1,nnp) 50 format(2x,2e22.10) c Since we have used MF=21, we have the eigenvalues c Check if we found purely imaginary or unstable eigenvalues write(6,51)IWORK(1),IWORK(2) 51 format(/ /,2x,'IWORK(1) : ',i4,2x,' IWORK(2) : ',i4) c Print eigenvalue information. N1=3+NNP+2*NNP*NNP+2*NNQ*NNQ N2=N1+NNP N3=N2+NNP+NNQ N4=N3+NNQ c the value of LWM is LWM = (MAXORD+1)*NEQN+20 if(MF.eq.21.or.MF.eq.20)MAXRD=6 if(MF.eq.11.or.MF.eq.10)MAXRD=13 c check against having set a non-default max-order c if(IOPT.eq.1.and.IWORK(5).ne.0)MAXRD=IWORK(5)+1 LWM=MAXRD*NEQN+20 c print the value of t where eigenvalues have been last computed write(6,52)TCJ 52 format(/ /,2x,'Last eigenvalues computation at TCJ: ',e22.10) c the closed loop eigenvalues are in A22-X*A12 for the c optimal control problems, ordered in descending magnitude write(6,53) 53 format(/ /,2x,' eigenvalues of A11(t)+A12(t)*X(t) '/' ') write(6,50)(RWORK(LWM+N1+I-1),RWORK(LWM+N2+I-1),I=1,NNP) c print final statistics on the run write(6,55) 55 format(/ /,2x,'Final statistics on this run '/' ') write(6,60)IWORK(11),IWORK(12),IWORK(13),RWORK(11),RWORK(12) 60 format 1 (2x,'nst= ',i5,' nfe=',i5,' nje=',i5,' hu=',e9.3,' hcur=',e9.3) stop end c subroutine PROBL(T,A,RARR,IARR) implicit real(a-h,o-z) dimension A(12,12) nn=12 do 10 j=1,nn do 10 i=1,nn A(i,j)=0.0e0 10 continue A(1,1)=20. A(4,1)=.744 A(5,1)=-.337 A(6,1)=-.02 A(1,7)=400. A(2,2)=25. A(4,2)=.032 A(5,2)=1.12 A(2,8)=625. A(6,3)=-.0386 A(4,4)=.154 A(5,4)=-.249 A(6,4)=.996 A(4,5)=.0042 A(5,5)=1. A(6,5)=.000295 A(4,6)=-1.54 A(5,6)=5.2 A(6,6)=.117 A(7,7)=-20. A(8,8)=-25. A(10,4)=1. A(7,10)=-.744 A(8,10)=-.032 A(10,10)=-.154 A(11,10)=-.0042 A(12,10)=1.54 A(11,5)=1. A(7,11)=.337 A(8,11)=-1.12 A(10,11)=.249 A(11,11)=-1. A(12,11)=-5.2 A(12,6)=1. A(7,12)=.02 A(9,12)=.0386 A(10,12)=-.996 A(11,12)=-.000295 A(12,12)=-.117 return end c c running this program produces the following output c on a Sun Sparc station with machine precision about 1.e-7 c c Current T and current solution are c c 0.1000000000E+11 c 0.3504407592E-02 0.4476269823E-03 c 0.3835782409E-02 -0.9914590418E-01 c -0.6655551028E-02 0.5426574498E-01 c 0.4476269823E-03 0.7874392904E-03 c 0.7484559901E-03 -0.2063505910E-01 c -0.1719598658E-01 0.4033283144E-01 c 0.3835782409E-02 0.7484559901E-03 c 0.1496508000E+08 -0.1110405475E+00 c -0.1308607869E-01 0.1001543552E+00 c -0.9914590418E-01 -0.2063505910E-01 c -0.1110405475E+00 0.2918051958E+01 c 0.3652955592E+00 -0.1965329170E+01 c -0.6655551028E-02 -0.1719598658E-01 c -0.1308607869E-01 0.3652955592E+00 c 0.3974425197E+00 -0.8232805133E+00 c 0.5426574498E-01 0.4033283144E-01 c 0.1001543552E+00 -0.1965329170E+01 c -0.8232805133E+00 0.5616596699E+01 c c c IWORK(1) : 0 IWORK(2) : 0 c c c eigenvalues of A11(t)+A12(t)*X(t) c c 0.3498644219E-06 0.0000000000E+00 c 0.7797037363E+00 0.1421359420E+01 c 0.7797037363E+00 -0.1421359420E+01 c 0.1647175312E+01 0.0000000000E+00 c 0.1998354912E+02 0.0000000000E+00 c 0.2497479439E+02 0.0000000000E+00 c c c Final statistics on this run c c nst= 33 nfe= 41 nje= 16 hu=0.435E+11 hcur=0.435E+11 c c -- driver 2 ended --