Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-01 - 43,50110/difeq.eng
There are 2 other files named difeq.eng in the archive. Click here to see a list.
3000'  NAME--DIFEQ
3010'
3020'  DESCRIPTION--SOLVES A SET OF N FIRST-ORDER DEFFERENTIAL EQUATIONS
3030'  BY THE FOURTH-ORDER RUNGE-KUTTA TECHNIQUE
3040'
3050'  SOURCE--A.O. CONVERSE, THAYER SCHOOL OF ENGINEERING
3060'
3070'  INSTRUCTIONS--THE FOLLOWING SYMBOLS ARE USED: X(I),T,G(I),N,H,
3080'  T1,H1,N1,N2,N3,N4,L(J,I),H2
3090'  
3100'  DX(I)/DT=G(I) FOR I=1,2,.....,N
3110'  WHERE THE G(I) ARE FUNCTIONS OF THE X(I) AND T.
3120'  N= THE NUMBER OF DEPENDENT VARIABLES; IT MUST BE <51
3130'  H= THE INCREMENT OF THE DEPENDENT VARIABLE
3140'  THE FUNCTIONS G(I) ARE TO BE INSERTED AFTER LINE 4000
3150'  AND BEFORE 4050 IN THE FOLLOWING FORM:
3160'  4000 LET G(1)= WHATEVER IT IS
3170'  4001 LET G(2)= WHATEVER IT IS
3175'    ETC.
3180'
3190'
3200'  *  *  *  *  *  *  *  MAIN PROGRAM  *  *  *  *  *  *  *  *  *  *  *
3210'
3220 DIM X(50)
3230 PRINT "HAVE YOU INSERTED THE 'LET G(I) = ...'"
3240 PRINT" STATEMENTS, YES OR NO"
3250 INPUTA1$
3260 IF A1$ = "YES" THEN 3290
3270PRINT "YOU MUST DO SO. LIST PROGRAM FOR EXPLANATION"
3280 GO TO 4100
3290PRINT"PLEASE TYPE THE NUMBER OF DEPENDENT VARIABLES FOLLOWED"
3300PRINT "BY A COMMA FOLLOWED BY THE INCREMENT SIZE OF THE INDEP."
3310 PRINT"VARIABLE FOLLOWED BY A COMMA FOLLOWED BY THE VALUE"
3320 PRINT"OF THE INDEP VAR. WHEN THE INTEGRATION IS TO STOP"
3330 INPUT N,H,T1
3340 PRINT "LIST THE INITIAL VALUES OF X(1),X(2),...,X(N) IN ORDER"
3350 MAT INPUT X
3360 PRINT "THE PROGRAM WILL PRINT T,X(N1),X(N2),X(N3), AND X(N4)"
3370 PRINT"EVERY H1 INCREMENTS. ENTER VALUES FOR H1,N1,N2,N3, AND"
3380 PRINT"N4 IN THAT ORDER AND SEPARATED BY COMMAS. ENTER VALUES FOR"
3390 PRINT" ALL TERMS EVEN IF N < 4. THE UNNEEDED VALUES WILL NOT"
3400 PRINT" AFFECT THE PROGRAM"
3410 INPUT H1,N1,N2,N3,N4
3420 PRINT
3430 PRINT
3440 PRINT"T","X(";N1;")","X(";N2;")","X(";N3;")","X(";N4;")"
3450 REM INITIAL VALUE OF THE INDEPENDENT VARIABLE IS 0.
3460 GO SUB 3820
3470 FOR I = 1 TO N
3480 LET L(1,I) = H*G(I)
3490 LET X(I) = X(I) + L(1,I)/2
3500 NEXT I
3510 LET T = T + H/2
3520 GO SUB 3820
3530 FOR I = 1 TO N
3540 LET L(2,I) = H*G(I)
3550 LET X(I) = X(I) -L(1,I)/2 + L(2,I)/2
3560 NEXT I
3570 GO SUB 3820
3580 FOR I = 1 TO N
3590 LET L(3,I) = H*G(I)
3600 LET X(I) = X(I) -L(2,I)/2 + L(3,I)
3610 NEXT I
3620 LET T = T + H/2
3630 GO SUB 3820
3640 FOR I = 1 TO N
3650 LET X(I) = X(I) - L(3,I) +(L(1,I)+2*L(2,I) +2*L(3,I) +H*G(I))/6
3660 NEXT I
3670 LET H2 = H2 + 1
3680 IF H2 < H1-.5 THEN 3780
3690 LET H2 = 0
3700 IF N < 1.5 THEN 3790
3710 IF N < 2.5 THEN 3770
3720 IF N < 3.5 THEN 3750
3730 PRINT T,X(N1),X(N2),X(N3),X(N4)
3740 GO TO 3800
3750 PRINT T, X(N1), X(N2), X(N3)
3760 GO TO 3800
3770 PRINT T,X(N1), X(N2)
3780 GO TO 3800
3790 PRINT T, X(N1)
3800 IF T > T1 THEN 4100 ' END OF INTEGRATION
3810 GO TO 3460
3820 REM THE G(I)'S ARE EVALUATED IN THIS SUBROUTINE
4000 LET G(1) = 1
4050 RETURN
4100 END