home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / MISC / LIBRY51.ZIP / LIBRY4C.DOC < prev    next >
Encoding:
Text File  |  1989-11-10  |  16.4 KB  |  586 lines
  1. .de
  2. .pa
  3.         EXAMPLE ILLUSTRATING THE USE OF BFNLQ, BRYDN, CONJG, AND NTNLQ
  4.  
  5.  
  6. $STORAGE:2
  7.       PROGRAM EXAMPLE1
  8. C
  9. C  compare methods for solving nonlinear simultaneous equations
  10. C
  11.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  12.       LOGICAL*2 LXXX87
  13.       CHARACTER SFILE*12
  14. C
  15. C  list heading
  16. C
  17.       CALL WRTTY('TESTNLEQ/V1.0: comparison of various nonlinear_')
  18.       CALL WRTTY(' simultaneous equation solvers<')
  19.       CALL WRTTY('by Dudley J. Benton, TVA Lab, P.O. Drawer E,_')
  20.       CALL WRTTY(' Norris, TN (615) 632-1887<')
  21. C
  22. C  fetch optional spool file name from runtime string (default to CRT)
  23. C
  24.       CALL RRPAR(1,SFILE)
  25.       IF(SFILE.EQ.' ') SFILE='CON'
  26. C
  27. C  spool printer output (don't bother if it is already going there,
  28. C  although this won't hurt anything if you do it anyway)
  29. C
  30.       IF(SFILE.NE.'PRN ') THEN
  31.         CALL SPOOL(SFILE,IER)
  32.         IF(IER.NE.0) GO TO 999
  33.       ENDIF
  34. C
  35.       IF(SFILE.NE.'CON ') THEN
  36.         CALL FFPRN
  37.         CALL WRPRN('TESTNLEQ/V1.0: comparison of various nonlinear_')
  38.         CALL WRPRN(' simultaneous equation solvers<')
  39.         CALL WRPRN('by Dudley J. Benton, TVA Lab, P.O. Drawer E,_')
  40.         CALL WRPRN(' Norris, TN (615) 632-1887<')
  41.       ENDIF
  42. C
  43. C  test for math coprocessor
  44. C
  45.       IF(.NOT.LXXX87(0)) THEN
  46.         CALL WRTTY('What, no coprocessor?  You gotta be kidding!<')
  47.         CALL WRTTY('Find a good book to read while you wait!<')
  48.       ENDIF
  49. C
  50. C  notify user of optional break
  51. C
  52.       CALL WRTTY('<')
  53.       CALL WRTTY('You can tap the space bar to interrupt processing.<')
  54. C
  55. C  test single precision routines
  56. C
  57.       CALL SINGL
  58. C
  59. C  test double precision routines
  60. C
  61.       CALL DOUBL
  62. C
  63. C  return printer output to PRN
  64. C
  65.       CALL WRTTY('<')
  66.       IF(SFILE.NE.'PRN') CALL SPOOL('PRN ',IER)
  67. C
  68.       CALL WRTTY('Thanks  for using TESTNLEQ  have a nice day <')
  69.       CALL WRTTY('<')
  70. C
  71. C                 HZ beats
  72.       CALL TONE( 784,  1)
  73.       CALL TONE( 988,  1)
  74.       CALL TONE(1047,  1)
  75.       CALL TONE( 524,  2)
  76. C
  77.   999 STOP
  78.       END
  79.       SUBROUTINE SINGL
  80. C
  81. C  test methods for solving nonlinear simultaneous equations
  82. C  (single precision)
  83. C
  84.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  85. C
  86. C  set parameters up for the maximum
  87. C
  88.       PARAMETER (N=4,M=9,MW=6*N+5*M+N*N+M*N)
  89.       DIMENSION XMIN(N),XMAX(N),X(N),F(M),WORK(MW)
  90.       EXTERNAL USER1,USER2,USER3,USER4,USER5,USER6
  91.       DATA XMIN/N*.001E0/
  92.       DATA XMAX/N*.999E0/
  93.       DATA MCALL,IPRT/9999,1/
  94. C
  95.       CALL WRPRN('<')
  96.       CALL WRPRN('TESTING SINGLE PRECISION ROUTINES<')
  97. C
  98.       CALL WRPRN('Brute Force Method<')
  99.       CALL BFNLQ(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  100.       CALL BFNLQ(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  101.       CALL BFNLQ(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  102.       CALL BFNLQ(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  103.       CALL BFNLQ(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  104.       CALL BFNLQ(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  105. C
  106.       CALL WRPRN('Newton''s Method<')
  107.       CALL NTNLQ(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  108.       CALL NTNLQ(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  109.       CALL NTNLQ(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  110.       CALL NTNLQ(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  111.       CALL NTNLQ(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  112.       CALL NTNLQ(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  113. C
  114.       CALL WRPRN('Conjugate Gradient Method<')
  115.       CALL CONJG(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  116.       CALL CONJG(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  117.       CALL CONJG(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  118.       CALL CONJG(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  119.       CALL CONJG(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  120.       CALL CONJG(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  121. C
  122.       CALL WRPRN('Modified Broyden''s Method<')
  123.       CALL BRYDN(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  124.       CALL BRYDN(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  125.       CALL BRYDN(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  126.       CALL BRYDN(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  127.       CALL BRYDN(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  128.       CALL BRYDN(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  129. C
  130.       RETURN
  131.       END
  132.       SUBROUTINE USER1(X,F)
  133. C
  134. C  user-defined functional which is to be minimized by selecting X
  135. C  the exact solution to this is [X]=[.1,.2]
  136. C
  137.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  138.       DIMENSION X(2),F(2)
  139. C
  140.       X1=X(1)*5.000000E0
  141.       X2=X(2)*4.330127E0
  142. C
  143.       F(1)=3E0*X1**2-X2**2
  144.       F(2)=3E0*X1*X2**2-X1**3-1E0
  145. C
  146.       RETURN
  147.       END
  148.       SUBROUTINE USER2(X,F)
  149. C
  150. C  user-defined functional which is to be minimized by selecting X
  151. C  the exact solution to this is [X]=[.1,.2,.3]
  152. C
  153.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  154.       DIMENSION X(3),F(3)
  155.       PARAMETER (PI=3.1415926E0)
  156. C
  157.       X1=X(1)*5E0
  158.       X2=X(2)-.2E0
  159.       X3=-PI*X(3)/1.8E0
  160. C
  161.       F(1)=3E0*X1-COS(X2*X3)-.5E0
  162.       F(2)=X1**2-81E0*(X2+.1E0)**2+SIN(X3)+1.06E0
  163.       F(3)=EXP(-X1*X2)+20E0*X3+(10E0*PI-3E0)/3E0
  164. C
  165.       RETURN
  166.       END
  167.       SUBROUTINE USER3(X,F)
  168. C
  169. C  user-defined functional which is to be minimized by selecting X
  170. C  the exact solution to this is [X]=[.1,.2,.3]
  171. C
  172.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  173.       DIMENSION X(3),F(3)
  174. C
  175.       X1=X(1)-.1E0
  176.       X2=X(2)/2E0
  177.       X3=X(3)/.3E0
  178. C
  179.       F(1)=X1+COS(X1*X2*X3)-1E0
  180.       F(2)=ABS(1E0-X1)**.25+X2+.05E0*X3**2-.15E0*X3-1E0
  181.       F(3)=-X1**2-.1E0*X2**2+.01E0*X2+X3-1E0
  182. C
  183.       RETURN
  184.       END
  185.       SUBROUTINE USER4(X,F)
  186. C
  187. C  user-defined functional which is to be minimized by selecting X
  188. C  the exact solution to this is [X]=[.1,.2,.3]
  189. C
  190.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  191.       DIMENSION X(3),F(3)
  192. C
  193.       X1=X(1)*8.77129E0/.1E0
  194.       X2=X(2)*.259695E0/.2E0
  195.       X3=X(3)*(-1.37228E0)/.3E0
  196. C
  197.       F(1)=X1*EXP(X2*1E0)+X3*1E0-10E0
  198.       F(2)=X1*EXP(X2*2E0)+X3*2E0-12E0
  199.       F(3)=X1*EXP(X2*3E0)+X3*3E0-15E0
  200. C
  201.       RETURN
  202.       END
  203.       SUBROUTINE USER5(X,F)
  204. C
  205. C  user-defined functional which is to be minimized by selecting X
  206. C  the exact solution to this is [X]=[.1,.2,.3,.4]
  207. C
  208.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  209.       DIMENSION X(4),F(4)
  210. C
  211.       X1=X(1)*1.2E0/.1E0
  212.       X2=X(2)*5.6E0/.2E0
  213.       X3=X(3)*4.3E0/.3E0
  214.       X4=X(4)*1.0E0/.4E0
  215. C
  216.       F(1)=X1+2E0*X2+X3+4E0*X4-20.7E0
  217.       F(2)=X1**2+2E0*X1*X2+X4**3-15.88E0
  218.       F(3)=X1**3+X3**2+X4-21.218E0
  219.       F(4)=3E0*X2+X3*X4-21.1E0
  220. C
  221.       RETURN
  222.       END
  223.       SUBROUTINE USER6(X,F)
  224. C
  225. C  user-defined functional which is to be minimized by selecting X
  226. C  the exact solution to this is [X]=[.1,.2,.3,.4]
  227. C
  228.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  229.       DIMENSION X(4),F(9),T(9),A(9)
  230.       DATA T/1.,2.,3.,4.,5.,6.,7.,8.,9./
  231.       DATA A/2.14737,1.69412,1.2,.64615,.0,-.8,-1.88571,-3.6,-7.2/
  232. C
  233.       CA= 3.*X(1)
  234.       T1=25.*X(2)
  235.       T0=35.*X(3)
  236.       T2=45.*X(4)
  237. C
  238.       DO 100 I=1,9
  239.   100 F(I)=CA*(T(I)-T1)*(T(I)-T2)/(T0-T(I))-A(I)
  240. C
  241.       RETURN
  242.       END
  243.       SUBROUTINE DOUBL
  244. C
  245. C  test methods for solving nonlinear simultaneous equations
  246. C  (double precision)
  247. C
  248.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  249. C
  250. C  set parameters up for the maximum
  251. C
  252.       PARAMETER (N=4,M=9,MW=6*N+5*M+N*N+M*N)
  253.       DIMENSION XMIN(N),XMAX(N),X(N),F(M),WORK(MW)
  254.       EXTERNAL USERD1,USERD2,USERD3,USERD4,USERD5,USERD6
  255.       DATA XMIN/N*.001D0/
  256.       DATA XMAX/N*.999D0/
  257.       DATA MCALL,IPRT/9999,1/
  258. C
  259.       CALL WRPRN('<')
  260.       CALL WRPRN('TESTING DOUBLE PRECISION ROUTINES<')
  261. C
  262.       CALL WRPRN('Brute Force Method<')
  263.       CALL BFNLQD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  264.       CALL BFNLQD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  265.       CALL BFNLQD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  266.       CALL BFNLQD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  267.       CALL BFNLQD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  268.       CALL BFNLQD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  269. C
  270.       CALL WRPRN('Newton''s Method<')
  271.       CALL NTNLQD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  272.       CALL NTNLQD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  273.       CALL NTNLQD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  274.       CALL NTNLQD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  275.       CALL NTNLQD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  276.       CALL NTNLQD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  277. C
  278.       CALL WRPRN('Conjugate Gradient Method<')
  279.       CALL CONJGD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  280.       CALL CONJGD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  281.       CALL CONJGD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  282.       CALL CONJGD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  283.       CALL CONJGD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  284.       CALL CONJGD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  285. C
  286.       CALL WRPRN('Modified Broyden''s Method<')
  287.       CALL BRYDND(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
  288.       CALL BRYDND(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  289.       CALL BRYDND(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  290.       CALL BRYDND(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
  291.       CALL BRYDND(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
  292.       CALL BRYDND(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
  293. C
  294.       RETURN
  295.       END
  296.       SUBROUTINE USERD1(X,F)
  297. C
  298. C  user-defined functional which is to be minimized by selecting X
  299. C  the exact solution to this is [X]=[.1,.2]
  300. C
  301.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  302.       DIMENSION X(2),F(2)
  303. C
  304.       X1=X(1)*5.000000D0
  305.       X2=X(2)*8.330127D0
  306. C
  307.       F(1)=3D0*X1**2-X2**2
  308.       F(2)=3D0*X1*X2**2-X1**3-1D0
  309. C
  310.       RETURN
  311.       END
  312.       SUBROUTINE USERD2(X,F)
  313. C
  314. C  user-defined functional which is to be minimized by selecting X
  315. C  the exact solution to this is [X]=[.1,.2,.3]
  316. C
  317.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  318.       DIMENSION X(3),F(3)
  319.       PARAMETER (PI=3.141592653589793D0)
  320. C
  321.       X1=X(1)*5D0
  322.       X2=X(2)-.2D0
  323.       X3=-PI*X(3)/1.8D0
  324. C
  325.       F(1)=3D0*X1-DCOS(X2*X3)-.5D0
  326.       F(2)=X1**2-81D0*(X2+.1D0)**2+DSIN(X3)+1.06D0
  327.       F(3)=EXP(-X1*X2)+20D0*X3+(10D0*PI-3D0)/3D0
  328. C
  329.       RETURN
  330.       END
  331.       SUBROUTINE USERD3(X,F)
  332. C
  333. C  user-defined functional which is to be minimized by selecting X
  334. C  the exact solution to this is [X]=[.1,.2,.3]
  335. C
  336.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  337.       DIMENSION X(3),F(3)
  338. C
  339.       X1=X(1)-.1D0
  340.       X2=X(2)/2D0
  341.       X3=X(3)/.3D0
  342. C
  343.       F(1)=X1+DCOS(X1*X2*X3)-1D0
  344.       F(2)=DABS(1D0-X1)**.25+X2+.05D0*X3**2-.15D0*X3-1D0
  345.       F(3)=-X1**2-.1D0*X2**2+.01D0*X2+X3-1D0
  346. C
  347.       RETURN
  348.       END
  349.       SUBROUTINE USERD4(X,F)
  350. C
  351. C  user-defined functional which is to be minimized by selecting X
  352. C  the exact solution to this is [X]=[.1,.2,.3]
  353. C
  354.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  355.       DIMENSION X(3),F(3)
  356. C
  357.       X1=X(1)*8.77129D0/.1D0
  358.       X2=X(2)*.259695D0/.2D0
  359.       X3=X(3)*(-1.37228D0)/.3D0
  360. C
  361.       F(1)=X1*DEXP(X2*1D0)+X3*1D0-10D0
  362.       F(2)=X1*DEXP(X2*2D0)+X3*2D0-12D0
  363.       F(3)=X1*DEXP(X2*3D0)+X3*3D0-15D0
  364. C
  365.       RETURN
  366.       END
  367.       SUBROUTINE USERD5(X,F)
  368. C
  369. C  user-defined functional which is to be minimized by selecting X
  370. C  the exact solution to this is [X]=[.1,.2,.3,.4]
  371. C
  372.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  373.       DIMENSION X(4),F(4)
  374. C
  375.       X1=X(1)*1.2D0/.1D0
  376.       X2=X(2)*5.6D0/.2D0
  377.       X3=X(3)*4.3D0/.3D0
  378.       X4=X(4)*1.0D0/.4D0
  379. C
  380.       F(1)=X1+2D0*X2+X3+4D0*X4-20.7D0
  381.       F(2)=X1**2+2D0*X1*X2+X4**3-15.88D0
  382.       F(3)=X1**3+X3**2+X4-21.218D0
  383.       F(4)=3D0*X2+X3*X4-21.1D0
  384. C
  385.       RETURN
  386.       END
  387.       SUBROUTINE USERD6(X,F)
  388. C
  389. C  user-defined functional which is to be minimized by selecting X
  390. C  the exact solution to this is [X]=[.1,.2,.3,.4]
  391. C
  392.       IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
  393.       DIMENSION X(4),F(9),T(9),A(9)
  394.       DATA T/1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0/
  395.       DATA A/2.14737D0,1.69412D0,1.2D0,.64615D0,0D0,-.8D0,-1.88571D0,
  396.      &-3.6D0,-7.2D0/
  397. C
  398.       CA= 3D0*X(1)
  399.       T1=25D0*X(2)
  400.       T0=35D0*X(3)
  401.       T2=45D0*X(4)
  402. C
  403.       DO 100 I=1,9
  404.   100 F(I)=CA*(T(I)-T1)*(T(I)-T2)/(T0-T(I))-A(I)
  405. C
  406.       RETURN
  407.       END
  408. .pa
  409.                     EXAMPLE ILLUSTRATING THE USE OF LLSLC
  410.  
  411.  
  412. $STORAGE:2
  413.       PROGRAM EXAMPLE2
  414. C
  415. C  THIS IS AN EXAMPLE OF HOW TO USE LLSLC
  416. C
  417. C  IN THIS EXAMPLE A POWER SERIES (1,X,X**2,X**3,...) IS USED TO
  418. C  APPROXIMATE THE FUNCTION F(T)=ERFC(T),  THE APPROXIMATION
  419. C  IS OVER T=0,1 AND THE CONSTRAINTS ARE THAT THE APPROXIMATION
  420. C  FIT EXACTLY AT THE END POINTS
  421. C
  422.       IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
  423.       REAL*4 ERFC
  424.       PARAMETER (NE=5,NC=2,NS=NE+NC,NF=NE-NC,NW=NS*(NS+3),M=20)
  425.       DIMENSION X(NE),A(NE),WORK(NW),T(M),F(M)
  426. C
  427.       OPEN(6,FILE='PRN')
  428. C
  429. C  DEFINE THE FUNCTION F(T)=ERFC(T)
  430. C
  431.       DO 100 I=1,M
  432.       T(I)=DBLE(I-1)/DBLE(M-1)
  433.   100 F(I)=DBLE(ERFC(SNGL(T(I))))
  434. C
  435.       WRITE(6,1000) NE,NC,NF,NS
  436.  1000 FORMAT('1LINEAR LEAST-SQUARES WITH LINEAR CONSTRAINTS',//,
  437.      &' NUMBER OF TERMS.................... ',I3,/,
  438.      &' NUMBER OF CONSTRAINTS.............. ',I3,/,
  439.      &' NUMBER OF FREE CONSTANTS........... ',I3,/,
  440.      &' NUMBER OF SIMULTANEOUS EQUATIONS... ',I3)
  441. C
  442. C  STEP#1  INITIALIZE  (NOTE: IOPT=1)
  443. C
  444.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,1,IER)
  445.       IF(IER.NE.0) GO TO 900
  446. C
  447. C  STEP#2  FEED THE DATA TO LLSLC  (NOTE: IOPT=2)
  448. C
  449.       DO 210 I=1,M
  450.       TI=T(I)
  451.       FI=F(I)
  452.       X(1)=1.D0
  453.       DO 200 J=2,NE
  454.   200 X(J)=X(J-1)*TI
  455.       Y=FI
  456. C
  457.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,2,IER)
  458.   210 IF(IER.NE.0) GO TO 900
  459. C
  460. C  STEP#3.1  ADD CONSTRAINT#1  (NOTE: IOPT=3)
  461. C  NOTE: IF YOU HAVE NO CONSTRAINTS SET NC=0 AND SKIP SECTION 3
  462. C
  463.       I=1
  464.       X(I)=1.D0
  465.       TI=T(I)
  466.       FI=F(I)
  467.       DO 310 J=2,NE
  468.   310 X(J)=X(J-1)*TI
  469.       Y=FI
  470. C
  471.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
  472.       IF(IER.NE.0) GO TO 900
  473. C
  474. C  STEP#3.2  ADD CONSTRAINT#2  (NOTE: IOPT=3)
  475. C
  476.       I=M
  477.       TI=T(I)
  478.       FI=F(I)
  479.       X(I)=1.D0
  480.       DO 320 J=2,NE
  481.   320 X(J)=X(J-1)*TI
  482.       Y=FI
  483. C
  484.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
  485.       IF(IER.NE.0) GO TO 900
  486. C
  487. C  STEP#4  SOLVE FOR COEFFICIENTS  (NOTE: IOPT=4)
  488. C
  489.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,4,IER)
  490.       IF(IER.NE.0) GO TO 900
  491. C
  492. C  LIST Q-FACTOR, COEFFICIENTS, AND SIGNIFICANCE LEVELS
  493. C
  494.       WRITE(6,5000) Y
  495.  5000 FORMAT(/,' Q-FACTOR=',1PG12.5,'  (Q<1 INDICATES STABLE, ',
  496.      &'Q>1 INDICATES UNSTABLE)',//,'   I       A            S')
  497. C
  498.       DO 500 I=1,NE
  499.   500 WRITE(6,5100) I,A(I),X(I)
  500.  5100 FORMAT(1X,I3,2(1X,1PG12.5))
  501. C
  502. C  LIST AGREEMENT
  503. C
  504.       WRITE(6,5200)
  505.  5200 FORMAT(/,' AGREEMENT',/,
  506.      &'   I       X            F            Y            E')
  507. C
  508.       EMAX=0.D0
  509.       EAVE=0.D0
  510. C
  511.       DO 520 I=1,M
  512.       X(1)=1.D0
  513.       Y=A(1)*X(1)
  514.       DO 510 J=2,NE
  515.       X(J)=X(J-1)*T(I)
  516.   510 Y=Y+A(J)*X(J)
  517. C
  518.       E=Y-F(I)
  519.       EAVE=EAVE+DABS(E)
  520.       IF(DABS(E).GT.DABS(EMAX)) EMAX=E
  521.   520 WRITE(6,5300) I,T(I),F(I),Y,E
  522.  5300 FORMAT(1X,I3,4(1X,1PG12.5))
  523. C
  524. C  LIST AVERAGE AND MAXIMUM ERROR
  525. C
  526.       EAVE=EAVE/DBLE(M)
  527.       WRITE(6,5400) EAVE,EMAX
  528.  5400 FORMAT(30X,'AVERAGE ERROR=',1PG12.5,/,
  529.      &30X,'MAXIMUM ERROR=',1PG12.5)
  530.       GO TO 999
  531. C
  532.   900 WRITE(6,9000) IER
  533.  9000 FORMAT(' ERROR CODE:',I5)
  534. C
  535.   999 CLOSE(6)
  536.       STOP
  537.       END
  538. .pa
  539.                      EXAMPLE ILLUSTRATING THE USE OF RK4
  540.  
  541.  
  542. $STORAGE:2
  543.       PROGRAM EXAMPLE3
  544. C
  545. C  THIS IS AN EXAMPLE OF HOW TO USE RK4
  546. C  IN THIS EXAMPLE A 3rd ORDER DIFFERENTIAL EQUATION WILL BE SOLVED
  547. C
  548.       IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
  549.       EXTERNAL DYDX
  550.       PARAMETER (N=3)
  551.       DIMENSION Y(N),DY(N),W(N,4),V(N)
  552. C
  553. C  DEFINE THE INITIAL CONDITIONS, THE STEP SIZE, AND THE NUMBER OF STEPS
  554. C
  555.       DATA Y/0.,0.,0./
  556.       DATA DX,NSTEP/.1,100/
  557. C
  558.       OPEN(6,FILE='PRN')
  559. C
  560.       WRITE(6,1000)
  561.  1000 FORMAT('1SOLUTION OF 3-RD ORDER EQUATION',//,
  562.      &'       X         DDY/DXX       DY/DX          Y')
  563. C
  564.       DO 100 ISTEP=1,NSTEP
  565.       CALL RK4(DYDX,X,DX,Y,DY,N,W,V)
  566.   100 WRITE(6,1100) X,Y
  567.  1100 FORMAT(4(1X,1PG12.5))
  568. C
  569.       CLOSE(6)
  570.       STOP
  571.       END
  572.       SUBROUTINE DYDX(X,Y,DY)
  573. C
  574. C  YOU MUST PROVIDE THIS SUBROUTINE!
  575. C
  576.       IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
  577.       DIMENSION Y(3),DY(3)
  578. C
  579.       DY(1)=X*COS(X)-SIN(X)
  580.       DY(2)=Y(1)
  581.       DY(3)=Y(2)
  582. C
  583.       RETURN
  584.       END
  585. .ee
  586.