home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power-Programmierung
/
CD1.mdf
/
fortran
/
library
/
library
/
libry4c.doc
< prev
next >
Wrap
Text File
|
1989-11-10
|
17KB
|
586 lines
.de
.pa
EXAMPLE ILLUSTRATING THE USE OF BFNLQ, BRYDN, CONJG, AND NTNLQ
$STORAGE:2
PROGRAM EXAMPLE1
C
C compare methods for solving nonlinear simultaneous equations
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
LOGICAL*2 LXXX87
CHARACTER SFILE*12
C
C list heading
C
CALL WRTTY('TESTNLEQ/V1.0: comparison of various nonlinear_')
CALL WRTTY(' simultaneous equation solvers<')
CALL WRTTY('by Dudley J. Benton, TVA Lab, P.O. Drawer E,_')
CALL WRTTY(' Norris, TN (615) 632-1887<')
C
C fetch optional spool file name from runtime string (default to CRT)
C
CALL RRPAR(1,SFILE)
IF(SFILE.EQ.' ') SFILE='CON'
C
C spool printer output (don't bother if it is already going there,
C although this won't hurt anything if you do it anyway)
C
IF(SFILE.NE.'PRN ') THEN
CALL SPOOL(SFILE,IER)
IF(IER.NE.0) GO TO 999
ENDIF
C
IF(SFILE.NE.'CON ') THEN
CALL FFPRN
CALL WRPRN('TESTNLEQ/V1.0: comparison of various nonlinear_')
CALL WRPRN(' simultaneous equation solvers<')
CALL WRPRN('by Dudley J. Benton, TVA Lab, P.O. Drawer E,_')
CALL WRPRN(' Norris, TN (615) 632-1887<')
ENDIF
C
C test for math coprocessor
C
IF(.NOT.LXXX87(0)) THEN
CALL WRTTY('What, no coprocessor? You gotta be kidding!<')
CALL WRTTY('Find a good book to read while you wait!<')
ENDIF
C
C notify user of optional break
C
CALL WRTTY('<')
CALL WRTTY('You can tap the space bar to interrupt processing.<')
C
C test single precision routines
C
CALL SINGL
C
C test double precision routines
C
CALL DOUBL
C
C return printer output to PRN
C
CALL WRTTY('<')
IF(SFILE.NE.'PRN') CALL SPOOL('PRN ',IER)
C
CALL WRTTY('Thanks for using TESTNLEQ have a nice day <')
CALL WRTTY('<')
C
C HZ beats
CALL TONE( 784, 1)
CALL TONE( 988, 1)
CALL TONE(1047, 1)
CALL TONE( 524, 2)
C
999 STOP
END
SUBROUTINE SINGL
C
C test methods for solving nonlinear simultaneous equations
C (single precision)
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
C
C set parameters up for the maximum
C
PARAMETER (N=4,M=9,MW=6*N+5*M+N*N+M*N)
DIMENSION XMIN(N),XMAX(N),X(N),F(M),WORK(MW)
EXTERNAL USER1,USER2,USER3,USER4,USER5,USER6
DATA XMIN/N*.001E0/
DATA XMAX/N*.999E0/
DATA MCALL,IPRT/9999,1/
C
CALL WRPRN('<')
CALL WRPRN('TESTING SINGLE PRECISION ROUTINES<')
C
CALL WRPRN('Brute Force Method<')
CALL BFNLQ(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQ(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQ(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQ(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQ(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQ(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Newton''s Method<')
CALL NTNLQ(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQ(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQ(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQ(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQ(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQ(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Conjugate Gradient Method<')
CALL CONJG(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL CONJG(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJG(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJG(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJG(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL CONJG(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Modified Broyden''s Method<')
CALL BRYDN(USER1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL BRYDN(USER2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDN(USER3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDN(USER4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDN(USER5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL BRYDN(USER6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
RETURN
END
SUBROUTINE USER1(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(2),F(2)
C
X1=X(1)*5.000000E0
X2=X(2)*4.330127E0
C
F(1)=3E0*X1**2-X2**2
F(2)=3E0*X1*X2**2-X1**3-1E0
C
RETURN
END
SUBROUTINE USER2(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(3),F(3)
PARAMETER (PI=3.1415926E0)
C
X1=X(1)*5E0
X2=X(2)-.2E0
X3=-PI*X(3)/1.8E0
C
F(1)=3E0*X1-COS(X2*X3)-.5E0
F(2)=X1**2-81E0*(X2+.1E0)**2+SIN(X3)+1.06E0
F(3)=EXP(-X1*X2)+20E0*X3+(10E0*PI-3E0)/3E0
C
RETURN
END
SUBROUTINE USER3(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(3),F(3)
C
X1=X(1)-.1E0
X2=X(2)/2E0
X3=X(3)/.3E0
C
F(1)=X1+COS(X1*X2*X3)-1E0
F(2)=ABS(1E0-X1)**.25+X2+.05E0*X3**2-.15E0*X3-1E0
F(3)=-X1**2-.1E0*X2**2+.01E0*X2+X3-1E0
C
RETURN
END
SUBROUTINE USER4(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(3),F(3)
C
X1=X(1)*8.77129E0/.1E0
X2=X(2)*.259695E0/.2E0
X3=X(3)*(-1.37228E0)/.3E0
C
F(1)=X1*EXP(X2*1E0)+X3*1E0-10E0
F(2)=X1*EXP(X2*2E0)+X3*2E0-12E0
F(3)=X1*EXP(X2*3E0)+X3*3E0-15E0
C
RETURN
END
SUBROUTINE USER5(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3,.4]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(4),F(4)
C
X1=X(1)*1.2E0/.1E0
X2=X(2)*5.6E0/.2E0
X3=X(3)*4.3E0/.3E0
X4=X(4)*1.0E0/.4E0
C
F(1)=X1+2E0*X2+X3+4E0*X4-20.7E0
F(2)=X1**2+2E0*X1*X2+X4**3-15.88E0
F(3)=X1**3+X3**2+X4-21.218E0
F(4)=3E0*X2+X3*X4-21.1E0
C
RETURN
END
SUBROUTINE USER6(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3,.4]
C
IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
DIMENSION X(4),F(9),T(9),A(9)
DATA T/1.,2.,3.,4.,5.,6.,7.,8.,9./
DATA A/2.14737,1.69412,1.2,.64615,.0,-.8,-1.88571,-3.6,-7.2/
C
CA= 3.*X(1)
T1=25.*X(2)
T0=35.*X(3)
T2=45.*X(4)
C
DO 100 I=1,9
100 F(I)=CA*(T(I)-T1)*(T(I)-T2)/(T0-T(I))-A(I)
C
RETURN
END
SUBROUTINE DOUBL
C
C test methods for solving nonlinear simultaneous equations
C (double precision)
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
C
C set parameters up for the maximum
C
PARAMETER (N=4,M=9,MW=6*N+5*M+N*N+M*N)
DIMENSION XMIN(N),XMAX(N),X(N),F(M),WORK(MW)
EXTERNAL USERD1,USERD2,USERD3,USERD4,USERD5,USERD6
DATA XMIN/N*.001D0/
DATA XMAX/N*.999D0/
DATA MCALL,IPRT/9999,1/
C
CALL WRPRN('<')
CALL WRPRN('TESTING DOUBLE PRECISION ROUTINES<')
C
CALL WRPRN('Brute Force Method<')
CALL BFNLQD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL BFNLQD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Newton''s Method<')
CALL NTNLQD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL NTNLQD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Conjugate Gradient Method<')
CALL CONJGD(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL CONJGD(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJGD(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJGD(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL CONJGD(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL CONJGD(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
CALL WRPRN('Modified Broyden''s Method<')
CALL BRYDND(USERD1,XMIN,XMAX,X,F,2,2,WORK,MW,MCALL,IPRT,IER)
CALL BRYDND(USERD2,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDND(USERD3,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDND(USERD4,XMIN,XMAX,X,F,3,3,WORK,MW,MCALL,IPRT,IER)
CALL BRYDND(USERD5,XMIN,XMAX,X,F,4,4,WORK,MW,MCALL,IPRT,IER)
CALL BRYDND(USERD6,XMIN,XMAX,X,F,4,9,WORK,MW,MCALL,IPRT,IER)
C
RETURN
END
SUBROUTINE USERD1(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(2),F(2)
C
X1=X(1)*5.000000D0
X2=X(2)*8.330127D0
C
F(1)=3D0*X1**2-X2**2
F(2)=3D0*X1*X2**2-X1**3-1D0
C
RETURN
END
SUBROUTINE USERD2(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(3),F(3)
PARAMETER (PI=3.141592653589793D0)
C
X1=X(1)*5D0
X2=X(2)-.2D0
X3=-PI*X(3)/1.8D0
C
F(1)=3D0*X1-DCOS(X2*X3)-.5D0
F(2)=X1**2-81D0*(X2+.1D0)**2+DSIN(X3)+1.06D0
F(3)=EXP(-X1*X2)+20D0*X3+(10D0*PI-3D0)/3D0
C
RETURN
END
SUBROUTINE USERD3(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(3),F(3)
C
X1=X(1)-.1D0
X2=X(2)/2D0
X3=X(3)/.3D0
C
F(1)=X1+DCOS(X1*X2*X3)-1D0
F(2)=DABS(1D0-X1)**.25+X2+.05D0*X3**2-.15D0*X3-1D0
F(3)=-X1**2-.1D0*X2**2+.01D0*X2+X3-1D0
C
RETURN
END
SUBROUTINE USERD4(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(3),F(3)
C
X1=X(1)*8.77129D0/.1D0
X2=X(2)*.259695D0/.2D0
X3=X(3)*(-1.37228D0)/.3D0
C
F(1)=X1*DEXP(X2*1D0)+X3*1D0-10D0
F(2)=X1*DEXP(X2*2D0)+X3*2D0-12D0
F(3)=X1*DEXP(X2*3D0)+X3*3D0-15D0
C
RETURN
END
SUBROUTINE USERD5(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3,.4]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(4),F(4)
C
X1=X(1)*1.2D0/.1D0
X2=X(2)*5.6D0/.2D0
X3=X(3)*4.3D0/.3D0
X4=X(4)*1.0D0/.4D0
C
F(1)=X1+2D0*X2+X3+4D0*X4-20.7D0
F(2)=X1**2+2D0*X1*X2+X4**3-15.88D0
F(3)=X1**3+X3**2+X4-21.218D0
F(4)=3D0*X2+X3*X4-21.1D0
C
RETURN
END
SUBROUTINE USERD6(X,F)
C
C user-defined functional which is to be minimized by selecting X
C the exact solution to this is [X]=[.1,.2,.3,.4]
C
IMPLICIT INTEGER*2(I-N),REAL*8(A-H,O-Z)
DIMENSION X(4),F(9),T(9),A(9)
DATA T/1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0/
DATA A/2.14737D0,1.69412D0,1.2D0,.64615D0,0D0,-.8D0,-1.88571D0,
&-3.6D0,-7.2D0/
C
CA= 3D0*X(1)
T1=25D0*X(2)
T0=35D0*X(3)
T2=45D0*X(4)
C
DO 100 I=1,9
100 F(I)=CA*(T(I)-T1)*(T(I)-T2)/(T0-T(I))-A(I)
C
RETURN
END
.pa
EXAMPLE ILLUSTRATING THE USE OF LLSLC
$STORAGE:2
PROGRAM EXAMPLE2
C
C THIS IS AN EXAMPLE OF HOW TO USE LLSLC
C
C IN THIS EXAMPLE A POWER SERIES (1,X,X**2,X**3,...) IS USED TO
C APPROXIMATE THE FUNCTION F(T)=ERFC(T), THE APPROXIMATION
C IS OVER T=0,1 AND THE CONSTRAINTS ARE THAT THE APPROXIMATION
C FIT EXACTLY AT THE END POINTS
C
IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
REAL*4 ERFC
PARAMETER (NE=5,NC=2,NS=NE+NC,NF=NE-NC,NW=NS*(NS+3),M=20)
DIMENSION X(NE),A(NE),WORK(NW),T(M),F(M)
C
OPEN(6,FILE='PRN')
C
C DEFINE THE FUNCTION F(T)=ERFC(T)
C
DO 100 I=1,M
T(I)=DBLE(I-1)/DBLE(M-1)
100 F(I)=DBLE(ERFC(SNGL(T(I))))
C
WRITE(6,1000) NE,NC,NF,NS
1000 FORMAT('1LINEAR LEAST-SQUARES WITH LINEAR CONSTRAINTS',//,
&' NUMBER OF TERMS.................... ',I3,/,
&' NUMBER OF CONSTRAINTS.............. ',I3,/,
&' NUMBER OF FREE CONSTANTS........... ',I3,/,
&' NUMBER OF SIMULTANEOUS EQUATIONS... ',I3)
C
C STEP#1 INITIALIZE (NOTE: IOPT=1)
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,1,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#2 FEED THE DATA TO LLSLC (NOTE: IOPT=2)
C
DO 210 I=1,M
TI=T(I)
FI=F(I)
X(1)=1.D0
DO 200 J=2,NE
200 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,2,IER)
210 IF(IER.NE.0) GO TO 900
C
C STEP#3.1 ADD CONSTRAINT#1 (NOTE: IOPT=3)
C NOTE: IF YOU HAVE NO CONSTRAINTS SET NC=0 AND SKIP SECTION 3
C
I=1
X(I)=1.D0
TI=T(I)
FI=F(I)
DO 310 J=2,NE
310 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#3.2 ADD CONSTRAINT#2 (NOTE: IOPT=3)
C
I=M
TI=T(I)
FI=F(I)
X(I)=1.D0
DO 320 J=2,NE
320 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#4 SOLVE FOR COEFFICIENTS (NOTE: IOPT=4)
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,4,IER)
IF(IER.NE.0) GO TO 900
C
C LIST Q-FACTOR, COEFFICIENTS, AND SIGNIFICANCE LEVELS
C
WRITE(6,5000) Y
5000 FORMAT(/,' Q-FACTOR=',1PG12.5,' (Q<1 INDICATES STABLE, ',
&'Q>1 INDICATES UNSTABLE)',//,' I A S')
C
DO 500 I=1,NE
500 WRITE(6,5100) I,A(I),X(I)
5100 FORMAT(1X,I3,2(1X,1PG12.5))
C
C LIST AGREEMENT
C
WRITE(6,5200)
5200 FORMAT(/,' AGREEMENT',/,
&' I X F Y E')
C
EMAX=0.D0
EAVE=0.D0
C
DO 520 I=1,M
X(1)=1.D0
Y=A(1)*X(1)
DO 510 J=2,NE
X(J)=X(J-1)*T(I)
510 Y=Y+A(J)*X(J)
C
E=Y-F(I)
EAVE=EAVE+DABS(E)
IF(DABS(E).GT.DABS(EMAX)) EMAX=E
520 WRITE(6,5300) I,T(I),F(I),Y,E
5300 FORMAT(1X,I3,4(1X,1PG12.5))
C
C LIST AVERAGE AND MAXIMUM ERROR
C
EAVE=EAVE/DBLE(M)
WRITE(6,5400) EAVE,EMAX
5400 FORMAT(30X,'AVERAGE ERROR=',1PG12.5,/,
&30X,'MAXIMUM ERROR=',1PG12.5)
GO TO 999
C
900 WRITE(6,9000) IER
9000 FORMAT(' ERROR CODE:',I5)
C
999 CLOSE(6)
STOP
END
.pa
EXAMPLE ILLUSTRATING THE USE OF RK4
$STORAGE:2
PROGRAM EXAMPLE3
C
C THIS IS AN EXAMPLE OF HOW TO USE RK4
C IN THIS EXAMPLE A 3rd ORDER DIFFERENTIAL EQUATION WILL BE SOLVED
C
IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
EXTERNAL DYDX
PARAMETER (N=3)
DIMENSION Y(N),DY(N),W(N,4),V(N)
C
C DEFINE THE INITIAL CONDITIONS, THE STEP SIZE, AND THE NUMBER OF STEPS
C
DATA Y/0.,0.,0./
DATA DX,NSTEP/.1,100/
C
OPEN(6,FILE='PRN')
C
WRITE(6,1000)
1000 FORMAT('1SOLUTION OF 3-RD ORDER EQUATION',//,
&' X DDY/DXX DY/DX Y')
C
DO 100 ISTEP=1,NSTEP
CALL RK4(DYDX,X,DX,Y,DY,N,W,V)
100 WRITE(6,1100) X,Y
1100 FORMAT(4(1X,1PG12.5))
C
CLOSE(6)
STOP
END
SUBROUTINE DYDX(X,Y,DY)
C
C YOU MUST PROVIDE THIS SUBROUTINE!
C
IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
DIMENSION Y(3),DY(3)
C
DY(1)=X*COS(X)-SIN(X)
DY(2)=Y(1)
DY(3)=Y(2)
C
RETURN
END
.ee