# 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
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