Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib10-01 - 43,50110/difeq1.clc
There are 2 other files named difeq1.clc in the archive. Click here to see a list.
```100'  NAME--DE-1OR
110'
120'  DESCRIPTION--SOLVE THE INITIAL VALUE PROBLEM
130'     Y' = F(X,Y)
140'     Y(X0) = Y0
150'
160'  SOURCE--RUNGA KUTTA METHOD (2ND ORDER ACCURACY)
170'
180'  INSTRUCTIONS--THE INTEGRATION STEP IS H AND THE SOLUTION IS
190'  TABULATED ON THE INTERVAL X0<= X <= B IN STEPS OF SIZE L.
200'
210'     350 DEF FNF(X,Y) = ........
220'     360 DATA X0,Y0,B,H,L
230'
240'  FOR MODIFICATIONS OF THE PRINT ROUTINE, FOR EXAMPLE, TO
250'  PLOT THE SOLUTION RATHER THAN TABULATE ITS VALUES, THE
260'  FOLLOWING LINES SHOULD BE MODIFIED:
270'     450 PRINT "VALUE OF X","VALUE OF Y"
280'     460 PRINT
290'     470 PRINT X,Y
300'     700 PRINT Q,Y1
310'
320'
330'  *  *  *  *  *  *   MAIN PROGRAM   *  *  *  *  *  *  *  *  *
340'
350 DEF FNF(X,Y)=-X/Y
360 DATA 0,1,1,.01,.1
380 LET X=X0
390 LET E=1/L
400 IF (B-X0)*H>0 THEN 420
410 LET H=-H
420 IF (B-X0)*L>0 THEN 440
430 LET L=-L
440 REM
450 PRINT "VALUE OF X","VALUE OF Y"
460 PRINT
470 PRINT X,Y
480 GOTO 520
490 LET P=Y+H*FNF(X,Y)
500 LET Y=(Y+P)/2+.5*H*FNF(X+H,P)
510 LET X=X+H
520 IF (X+H-B)*SGN(L)+1E-7 >=0 THEN 550
530 LET R=X+H
540 GOTO 560
550 LET R=B
560 LET A=INT(E*R)/E + L*(SGN(H)-1)/2
570 IF (A+L-R)*SGN(L)-1E-7 >=0 THEN 590
580 LET A=A+L
590 IF (X-A)*SGN(L)+1E-7 >= 0 THEN 620
600 LET Q=A
610 GOSUB 680
620 IF R=B THEN 640
630 GOTO 490
640 IF ABS(R-A)<.5*1E-6 THEN 670
650 LET Q=R
660 GOSUB 680
670 STOP
680 LET P1=Y+(Q-X)*FNF(X,Y)
690 LET Y1=(Y+P1)/2+.5*(Q-X)*FNF(Q,P1)
700 PRINT Q,Y1
710 RETURN
720 END

```