OptDes via Quasilinearization

This programmed example was taken from the Appendix of the textbook, Quasilinearization and Nonlinear Boundary-Value Problems This is only one of many example problems from this text, which are solved by quasilinearization in a similar unified program structure, where the context is that of the algorithm, not a problem statement.

         PROGRAM OPTDES(INPUT,OUTPUT)
C DESIGN PROBLEM
      COMMON N1,KMAX,NP,NMAX,HGRID,CX,T,X
     ~       Y,U,V,H,PREV,HP,A,B,C,NEQ
      DIMENSION T(250),W(4,500),H(4,4,500)
     ~          P(4,500),A(2,2),B(2),C(4),PREV(4)
C       X, STATE VARIABLE
C       Y, CONTROL VARIABLE
C       U(=MU), LAGRANGE MULTIPLIER
C       V(=A), DESIGN PARAMETER
C
   1  CALL INPUT
   2  CALL START
C
C                    K ITERATIONS
C
   3  DO 19 K=1,KMAX
C                    INITIALIZE
      DO 5 I =1,250
   5  T(I)=0.
      T(2)=0.
      T(3)=HGRID
      T(4)=1.
      T(9)=1.
      T(14)=1.
      N=0
      PRINT 111,(T(L),L=4,43)
      NEC=20
      CALL INT(T,NEQ,N1,0.,0.,0.,0.,0.,0.)
C           INTEGRATE P'S AND H'S STORING GRID VALUES
      N=0
   8  DO 11 N=1,NP
      N=N+1
        X=W(1,N)
        Y=W(2,N)
        U=W(3,N)
        V=W(4,N)
        CALL INTM
        L=3
        DO 9 I=1,4
           DO 9 J=1,4
           L=L+1
   9    H(I,J,N)=T(L)
        DO 10 J=1,4
          L=L+1
  10    P(J,N)=T(L)
  11  CONTINUE
C
   6  PRINT 111,N,(T(L),L=4,43
      PRINT 111,N,(H(I,J,N),J=1,4),(P(J,N),J=1,4)
 111  FORMAT(1H019XI10,4E20.6/30X4E20.6)
      IF(N-NMAX)8,7,7
   7  PRINT 40,K
  40    FORMAT(1H0/65X9HITERATION,I3)
C
C                 EVALUATE CONSTANTS
C
      N=NMAX
      A(1,1)=H(2,2,N)
      A(2,1)=H(2,3,N)
      A(1,2)=H(3,2,N) + H(4,2,N)
      A(2,2)=H(3,3,N) + H(4,3,N)
C
      B(1)=-P(2,N) - CX*H(1,2,N)
      B(2)=-P(3,N) - CX*H(1,3,N)
C
      DET=A(1,1)*A(2,2) - A(2,1)*A(1,2)
      G=ABSF(DET)
      IF(G-.000001)200,200,12
C
 200  PRINT 201,DET,((A(I,J),J=1,2),B(I),I=1,2)
 201    FORMAT(1H020X12HDETERMINANT=,E16.6/(40X3E20.6))
      CALL EXIT
  12  C(1)=CX
      C(2)=(B(1)*A(2,2) - B(2)*A(1,2))/DET
      C(3)=(B(2)*A(1,1) - B(1)*A(2,1))/DET
      C(4)=C(3)
      GRID=0.
      PRINT 45, GRID,(C(I),I=1,4)
  45    FORMAT(1H026XHGRID,14X1HX,19X1HY,18X2HMU,19X1HA/
     *  20XF10.4,4E20.6)
C
C                 NEW VALUES OF VARIABLES
      N=0
  13  DO 14 M=1,NP
        N=N+1
        DO 14 I=1,4
        W(I,N)=P(I,N)
        DO 14 J=1,4
  14    W(I,N)=W(I,N) +C(J)*H(J,I,N)
      FN=N
      GRID=FN*HGRID
      PRINT 50, GRID,(W(I,N),I=1,4)
      IF(N-NMAX)13,15,15
C
C                 COMPARE C'S, PRESENT AND PREVIOUS
  15  DO 16 1,5
        G=ABSF(C(I)-PREV(I))
        IF(G-.000001)16,16,17
  16  CONTINUE
          GO TO 1
  17  DO 18 I=1,4
  18  PREV(I)=C(I)
C
  19  CONTINUE
          GO TO 1
C 50    FORMAT(20XF10.4,4E20.6)
      END
      SUBROUTINE INPUT
      COMMON N1,KMAX,HGRID,CX,T,X,Y,U,V,W,PREV,H,P,A,B,C
      DIMENSION T(250),W(4,500),H(4,4,500),P(4,500),A(2,2),B(2),C(4),
     * PREV(4)
C          N1, INTEGRATION OPTION
C          KMAX, NUMBER OF ITERATIONS
C          NP, NUMBER OF GRIDS BEFORE PRINT-OUT
C          NMAX, TOTAL NUMBER OF INTEGRATIONS IN INTERVAL
C          HGRID, GRID SIZE
C          CX, X(0)
C
      READ 110,N1,KMAX,NP,NMAX,HGRID,CX
      PRINT ID,N1,KMAX,NP,NMAX,HGRID,CX
 110   FORMAT(4I4,2E8.2)
  10   FORMAT(1H120X14HDESIGN PROBLEM//
     1 38X2HN1,6X4HKMAX,8X2HNP,6X4HNMAX,13X5HHGRID,16X4HX(0)/
     2 30X4I10,2E20.6)
        RETURN
      END

      SUBROUTINE START
       COMMON N1,KMAX,NP,NMAX,HGRID,CX,T,X,Y,U,V,W,PREV,H,P,A,B,C
       DIMENSION T(250),W(4,500),H(4,4,500),P(4,500),A(2,2),B(2),C(4),
     1  PREV(4)
C           INITIALLY T=0, X(T)=CX, Y(T)=0. A(T)=0
      DO 1 N=1,NMAX
      W(1,N)=CX
      W(2,N)=0.
      W(3,N)=0.
   1  W(4,N)=0.
C
      DO 2 I=1,4
   2  PREV(I)=W(I,1)
C
      K=0
      PRINT 40,K
      GRID=0.
      PRINT 45,GRID,(W(I,1),I=1,4)
  40    FORMAT(1H0/65X9HITERATION,I3)
  45    FORMAT(1H026X4HGRID,14X1HX,19X1HY,18X2HMU,19X1HA/20XF10.4,4E20.6
     1  /1H019X22HFOR ALL VALUES OF TIME)
          RETURN
      END

      SUBROUTINE DAUX
       COMMON N1,KMAX,NP,NMAX,HGRID,CX,T,X,Y,U,V,W,PREV,H,P,A,B,C,NEQ
       DIMENSION T(250),W(4,500),H(4,4,500),P(4,500),A(2,2),8(2),C(4),
     1  PREV(4)
C
      L=20
      DO 1 I=1,4
        L=L+4
      T(L)=-V*T(L-20) + T(L-19) -X*T(L-17)
      T(L+1)=T(L-20) + V*T(L-19) + Y(T(L-17)
      T(L+2)=Y*T(L-20) + X*T(L-19)
   1  T(L+3)=0.
C
      T(40)=-V*T(20) + T(21) - X*T(23) +V*X
      T(41)=T(20) + V*T(21) + Y*T(23) - V*Y
      T(42)=Y*T(20) + X*T(21) - X*Y
      T(43)=0.
        RETURN
      END