C PCPRB2.FOR 27 February 1991 C PROGRAM PCPRB2 C C The aircraft stability problem. C C The variables X(I) are control parameters for an aircraft under C the guidance of a pilot: C C C X(1) = Roll C X(2) = Pitch C X(3) = Yaw C X(4) = Angle of attack C X(5) = Sideslip C X(6) = Elevator C X(7) = Aileron C X(8) = Rudder C C C The function F is of the following form: C C C For I=1 to 5, C C F(I)=Sum(J=1,8) Sum(K=1,8) B(I,J)*X(J)+PHI(I,J,K)*X(J)*X(K) C C F(6)=X(IFIX1)-VAL1 C F(7)=X(IFIX2)-VAL2 C C C Options: C C C ICHOOZ=1 IFIX1=6, VAL1=-.050, Rheinboldt BARRAY. C ICHOOZ=2 IFIX1=6, VAL1=-.008, Rheinboldt BARRAY. C ICHOOZ=3 IFIX1=6, VAL1=0.000, Rheinboldt BARRAY. C ICHOOZ=4 IFIX1=6, VAL1=0.050, Rheinboldt BARRAY. C ICHOOZ=5 IFIX1=6, VAL1=0.100, Rheinboldt BARRAY. C ICHOOZ=6 IFIX1=6, VAL1=-0.01250, Rheinboldt BARRAY. (Bifurcation) C C ICHOOZ=11 IFIX1=6, VAL1=-.050, Melhem BARRAY. C ICHOOZ=12 IFIX1=6, VAL1=-.008, Melhem BARRAY. C ICHOOZ=13 IFIX1=6, VAL1=0.000, Melhem BARRAY. C ICHOOZ=14 IFIX1=6, VAL1=0.050, Melhem BARRAY. C ICHOOZ=15 IFIX1=6, VAL1=0.100, Melhem BARRAY. C ICHOOZ=16 IFIX1=6, VAL1=-0.01250, Melhem BARRAY. (Bifurcation) C C For all choices of ICHOOZ, IFIX2=8, VAL2=0.0. C If ICHOOZ=6 or 16, LIM=0, otherwise LIM=7. C C C Limit points: C C C Melhem lists the following limit points in X7. Note that Melhem C has BARRAY(4,1)=1.0, BARRAY(4,2)=0.0: C C C X1 X2 X3 X4 X5 X6 X7 X8 C C -2.9691 0.83074 -0.072748 0.41029 -0.26880 -.05 0.50919 0.0 C -2.8158 -0.17481 -0.089469 0.026319 0.070951 -.008 0.20442 0.0 C -3.7571 -0.64911 -0.39350 0.091810 0.19685 -.008 -.0038238 0.0 C -4.1637 0.092284 -0.092610 0.022402 -0.017106 -.008 0.37823 0.0 C -2.5839 -0.22128 -0.054079 0.013524 0.090871 0.0 0.18608 0.0 C -3.9007 -1.1421 -0.57863 0.13284 0.32685 0.0 -0.50703 0.0 C -2.3610 -0.72360 0.032739 -0.039108 0.29347 0.05 0.29272 0.0 C -2.2982 1.4033 0.063244 -0.079383 0.58336 0.10 0.58336 0.0 C C C Rheinboldt lists the following limit points, using C BARRAY(4,1)=0.0, BARRAY(4,2)=1.0: C C C X1 X2 X3 X4 X5 X6 X7 X8 C C 2.9648 0.82556 0.07366 0.041309 0.26734 -0.050 -.050481 0.0 C 2.8173 -0.17628 0.08992 0.026429 -0.07147 -0.008 -.204973 0.0 C 3.7579 -0.65541 0.38658 0.092520 -0.19867 -0.008 0.006200 0.0 C 4.1638 0.08913 0.09480 0.022888 0.16232 -0.008 -.377660 0.0 C 2.5873 -0.22354 0.05468 0.013676 -0.09168 0.000 -.186908 0.0 C 3.9005 -1.14815 0.58156 0.133516 -0.32858 0.000 .510158 0.0 C 2.3639 -0.72974 -0.31604 -0.038785 -0.29583 0.050 -.295772 0.0 C 2.2992 -1.41023 -0.06184 -0.079009 -0.58629 0.100 -.689717 0.0 C C C Bifurcation points: C C C Rheinboldt, using BARRAY(4,1)=0.0, BARRAY(4,2)=1.0, lists: C C C 4.482 0.1632 0.02373 0.006205 0.03527 -0.0006177 -0.3986 0.0 C 3.319 -0.1869 0.1605 0.04379 -0.06888 -0.01250 -0.2374 0.0 C 4.466 0.1467 0.04045 0.009777 0.03089 -0.006129 -0.3995 0.0 C -3.325 0.1880 -0.1614 0.04395 0.06911 -0.01247 0.2367 0.0 C C References: C C A A Schy, M E Hannah, C Prediction of Jump Pheonomena in Roll-Coupled Maneuvers of Airplanes, C Journal of Aircraft, 14, 1977, Pages 375-382. C C J W Young, A A Schy, K G Johnson C Prediction of Jump Phenomena in Aircraft Maneuvers, Including Nonlinear C Aerodynamic Effects, C Journal of Guidance and Control 1, 1978, Pages 26-31. C INTEGER LIW INTEGER LRW INTEGER NVAR C PARAMETER (NVAR=8) PARAMETER (LIW=29+NVAR) PARAMETER (LRW=29+NVAR*(NVAR+4)) C EXTERNAL DENSLV EXTERNAL FPAIR EXTERNAL FXAIR EXTERNAL PITCON C REAL BARRAY(5,8) REAL HMAX INTEGER I INTEGER ICHOOZ INTEGER IERROR INTEGER IFIX1 INTEGER IFIX2 INTEGER IPAR(2) INTEGER IWORK(LIW) INTEGER J INTEGER K INTEGER LIM INTEGER LOUNIT CHARACTER NAME*12 REAL RPAR(42) REAL RWORK(LRW) REAL VAL1 REAL VAL2 REAL XR(NVAR) C LOUNIT=6 C C Set options. C ICHOOZ=1 IFIX1=6 IFIX2=8 C C Initialize work arrays. C DO 4 I=1,LIW IWORK(I)=0 4 CONTINUE DO 5 I=1,LRW RWORK(I)=0.0 5 CONTINUE C C Set VAL1, VAL2 based on ICHOOZ C IF(ICHOOZ.EQ.1.OR.ICHOOZ.EQ.11)VAL1=-.050 IF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.12)VAL1=-.008 IF(ICHOOZ.EQ.3.OR.ICHOOZ.EQ.13)VAL1=0.000 IF(ICHOOZ.EQ.4.OR.ICHOOZ.EQ.14)VAL1=0.050 IF(ICHOOZ.EQ.5.OR.ICHOOZ.EQ.15)VAL1=0.100 IF(ICHOOZ.EQ.6.OR.ICHOOZ.EQ.16)VAL1=-0.01250 VAL2=0.0 C C Set some parameters C IPAR(1)=IFIX1 IPAR(2)=IFIX2 IF(ICHOOZ.NE.6.AND.ICHOOZ.NE.16)THEN HMAX=0.4 ELSE HMAX=0.2 ENDIF IF(ICHOOZ.NE.6.AND.ICHOOZ.NE.16)THEN LIM=7 ELSE LIM=0 ENDIF C C Set starting point C DO 10 I=1,NVAR XR(I)=0.0 10 CONTINUE XR(IFIX1)=VAL1 XR(IFIX2)=VAL2 C C Set BARRAY C BARRAY(1,1)=-3.933 BARRAY(2,1)=0.0 BARRAY(3,1)=0.002 IF(ICHOOZ.LT.10)THEN BARRAY(4,1)=0.0 ELSE BARRAY(4,1)=1.0 ENDIF BARRAY(5,1)=0.0 BARRAY(1,2)=0.107 BARRAY(2,2)=-0.987 BARRAY(3,2)=0.0 IF(ICHOOZ.LT.10)THEN BARRAY(4,2)=1.0 ELSE BARRAY(4,2)=0.0 ENDIF BARRAY(5,2)=0.0 BARRAY(1,3)=0.126 BARRAY(2,3)=0.0 BARRAY(3,3)=-0.235 BARRAY(4,3)=0.0 BARRAY(5,3)=-1.0 BARRAY(1,4)=0.0 BARRAY(2,4)=-22.95 BARRAY(3,4)=0.0 BARRAY(4,4)=-1.0 BARRAY(5,4)=0.0 BARRAY(1,5)=-9.99 BARRAY(2,5)=0.0 BARRAY(3,5)=5.67 BARRAY(4,5)=0.0 BARRAY(5,5)=-0.196 BARRAY(1,6)=0.0 BARRAY(2,6)=-28.37 BARRAY(3,6)=0.0 BARRAY(4,6)=-0.168 BARRAY(5,6)=0.0 BARRAY(1,7)=-45.83 BARRAY(2,7)=0.0 BARRAY(3,7)=-0.921 BARRAY(4,7)=0.0 BARRAY(5,7)=-0.0071 BARRAY(1,8)=-7.64 BARRAY(2,8)=0.0 BARRAY(3,8)=-6.51 BARRAY(4,8)=0.0 BARRAY(5,8)=0.0 C C Copy BARRAY, VAL1, VAL2 into RPAR. C K=0 DO 30 I=1,5 DO 29 J=1,NVAR K=K+1 RPAR(K)=BARRAY(I,J) 29 CONTINUE 30 CONTINUE K=K+1 RPAR(K)=VAL1 K=K+1 RPAR(K)=VAL2 C C Set work array input C C IWORK(1)=0 ; This is a startup C IWORK(2)=7 ; Use X(7) for initial parameter C IWORK(3)=0 ; Program may choose parameter index C IWORK(4)=0 ; Update jacobian every newton step C IWORK(5)=0 ; Seek no target values. C IWORK(6)=LIM ; Seek limit points in X(LIM) C IWORK(7)=1 ; Control amount of output. C IWORK(8)=6 ; Output unit C IWORK(9)=0 ; Jacobian choice. C IWORK(1)=0 IWORK(2)=7 IWORK(3)=0 IWORK(4)=0 IWORK(5)=0 IWORK(6)=LIM IWORK(7)=1 IWORK(8)=LOUNIT IWORK(9)=0 C C RWORK(1)=0.0001 ; Absolute error tolerance C RWORK(2)=0.0001 ; Relative error tolerance C RWORK(3)=0.0001 ; Minimum stepsize C RWORK(4)=HMAX ; Maximum stepsize C RWORK(5)=0.1 ; Starting stepsize C RWORK(6)=-1.0 ; Starting direction C RWORK(7)=0.0 ; Target value C RWORK(1)=0.0001 RWORK(2)=0.0001 RWORK(3)=0.0001 RWORK(4)=HMAX RWORK(5)=0.1 RWORK(6)=-1.0 RWORK(7)=0.0 C WRITE(6,*)' ' WRITE(6,*)'PCPRB2:' WRITE(6,*)'PITCON sample program.' WRITE(6,*)' ' WRITE(6,*)'The Aircraft Stability Problem' WRITE(6,*)'Using option ',ICHOOZ WRITE(6,*)'Fix variable ',IFIX1,' at ',VAL1 WRITE(6,*)'Fix variable ',IFIX2,' at ',VAL2 IF(ICHOOZ.LT.10)THEN WRITE(6,*)'Rheinboldt version of BARRAY' ELSE WRITE(6,*)'Melhem version of BARRAY' ENDIF WRITE(6,*)' ' WRITE(6,*)'Step Type of Point ' WRITE(6,*)' ' I=0 NAME='Start point ' WRITE(6,'(1X,I3,2X,A12,2X,8F7.3)')I,NAME,(XR(J),J=1,NVAR) C DO 40 I=1,50 CALL PITCON(FPAIR,RPAR,FXAIR,IERROR,IPAR,IWORK,LIW, * NVAR,RWORK,LRW,XR,DENSLV) C IF(IWORK(1).EQ.1)THEN NAME='Corrected ' ELSEIF(IWORK(1).EQ.2)THEN NAME='Continuation ' ELSEIF(IWORK(1).EQ.3)THEN NAME='Target point ' ELSEIF(IWORK(1).EQ.4)THEN NAME='Limit point ' ENDIF C WRITE(6,'(1X,I3,2X,A12,2X,8F7.3)')I,NAME,(XR(J),J=1,NVAR) C IF(IERROR.NE.0)THEN WRITE(6,*)'Iteration halted, IERROR=',IERROR GO TO 50 ENDIF C 40 CONTINUE 50 CONTINUE C WRITE(6,*)' ' WRITE(6,*)'Jacobians: ',IWORK(19) WRITE(6,*)'Factorizations:',IWORK(20) WRITE(6,*)'Solves: ',IWORK(21) WRITE(6,*)'Functions: ',IWORK(22) C STOP END SUBROUTINE FXAIR(NVAR,RPAR,IPAR,X,FX,IERROR) C C*********************************************************************** C INTEGER NVAR C REAL FX(NVAR) INTEGER I INTEGER IERROR INTEGER IFIX1 INTEGER IFIX2 INTEGER IPAR(2) INTEGER J INTEGER K REAL PHI REAL RPAR(42) REAL VAL1 REAL VAL2 REAL X(NVAR) C C Compute linear terms C K=0 DO 20 I=1,5 FX(I)=0.0 DO 10 J=1,8 K=K+1 FX(I)=FX(I)+RPAR(K)*X(J) 10 CONTINUE 20 CONTINUE C C Compute nonlinear terms C PHI=-.727*X(2)*X(3)+8.39*X(3)*X(4)-684.4*X(4)*X(5)+63.5*X(4)*X(7) FX(1)=FX(1)+PHI C PHI=.949*X(1)*X(3)+.173*X(1)*X(5) FX(2)=FX(2)+PHI C PHI=-.716*X(1)*X(2)-1.578*X(1)*X(4)+1.132*X(4)*X(7) FX(3)=FX(3)+PHI C PHI=-X(1)*X(5) FX(4)=FX(4)+PHI C PHI=X(1)*X(4) FX(5)=FX(5)+PHI C C Two function values restrict two variables: C IFIX1=IPAR(1) VAL1=RPAR(41) FX(6)=X(IFIX1)-VAL1 C IFIX2=IPAR(2) VAL2=RPAR(42) FX(7)=X(IFIX2)-VAL2 C RETURN END SUBROUTINE FPAIR(NVAR,RPAR,IPAR,X,FPRIME,IERROR) C C*********************************************************************** C INTEGER NVAR C REAL FPRIME(NVAR,NVAR) INTEGER I INTEGER IERROR INTEGER IFIX1 INTEGER IFIX2 INTEGER IPAR(2) INTEGER J INTEGER K REAL RPAR(42) REAL X(NVAR) C C Linear part of F C K=0 DO 20 I=1,5 DO 10 J=1,8 K=K+1 FPRIME(I,J)=RPAR(K) 10 CONTINUE 20 CONTINUE C DO 40 I=6,7 DO 30 J=1,8 FPRIME(I,J)=0.0 30 CONTINUE 40 CONTINUE C IFIX1=IPAR(1) IFIX2=IPAR(2) FPRIME(6,IFIX1)=1.0 FPRIME(7,IFIX2)=1.0 C C Nonlinear part of F C FPRIME(1,2)=FPRIME(1,2)-.727*X(3) FPRIME(1,3)=FPRIME(1,3)-.727*X(2)+8.39*X(4) FPRIME(1,4)=FPRIME(1,4)+8.39*X(3)-684.4*X(5)+63.5*X(7) FPRIME(1,5)=FPRIME(1,5)-684.4*X(4) FPRIME(1,7)=FPRIME(1,7)+63.5*X(4) C FPRIME(2,1)=FPRIME(2,1)+.949*X(3)+.173*X(5) FPRIME(2,3)=FPRIME(2,3)+.949*X(1) FPRIME(2,5)=FPRIME(2,5)+.173*X(1) C FPRIME(3,1)=FPRIME(3,1)-.716*X(2)-1.578*X(4) FPRIME(3,2)=FPRIME(3,2)-.716*X(1) FPRIME(3,4)=FPRIME(3,4)-1.578*X(1)+1.132*X(7) FPRIME(3,7)=FPRIME(3,7)+1.132*X(4) C FPRIME(4,1)=FPRIME(4,1)-X(5) FPRIME(4,5)=FPRIME(4,5)-X(1) C FPRIME(5,1)=FPRIME(5,1)+X(4) FPRIME(5,4)=FPRIME(5,4)+X(1) C RETURN END