145 lines
7.9 KiB
Plaintext
145 lines
7.9 KiB
Plaintext
[Demo of EDSAC library subroutine G1: Runge-Kutta solution of differential equations.
|
|
Full description is in Wilkes, Wheeler & Gill, 1951 edn, pages 32-34, 86-87, 132-134.
|
|
|
|
Before using G1, we need to fix n, m, a, b, c, d, as defined in WWG pages 86-87:
|
|
n = number of equations (2 for the Rosetta Code example).
|
|
2^m = multiplier for the hy', as large as possible without causing numeric overflow;
|
|
with the scaling chosen here, m = 5.
|
|
Variables y are stored in n consecutive long locations, the last of which is aD.
|
|
Scaled derivatives (2^m)hy' in n consecutive long locations, the last of which is bD.
|
|
G1 uses working variables in n consecutive long locations, the last of which is cD.
|
|
d = address of user-supplied auxiliary subroutine, which calculates the (2^m)hy'.
|
|
|
|
For convenience, keep G1 and its storage together. Start at (say) 400 and place:
|
|
variables y at 400D, 402D;
|
|
scaled derivatives at 404D, 406D;
|
|
workspace for G1 at 408D, 410D;
|
|
G1 itself at 412.
|
|
If the base address is placed in location 51 at load time, all the above
|
|
addresses can be accessed via the G parameter:]
|
|
T 51 K
|
|
P 400 F
|
|
[Now set up the 6 preset parameters specified in WWG:]
|
|
T 45 K
|
|
P 2#G [H parameter: P a D]
|
|
P 4 F [N parameter: P 2n F]
|
|
P 4 F [M parameter: P (b-a) F, or V (2048-a+b) F if a > b]
|
|
P 4 F [& parameter: P (c-b) F, or V (2048-b+c) F if b > c]
|
|
P 8 F [L parameter: P 2^(m-2) F]
|
|
P 300 F [X parameter: P d F]
|
|
[For other addresses in the program we can optionally use some more parameters:]
|
|
T 52 K
|
|
P 120 F [A parameter: main routine]
|
|
P 56 F [B parameter: print subroutine P1 from EDSAC library]
|
|
P 350 F [C parameter: constants for Rosetta code example]
|
|
P 78 F [V parameter: square root subroutine]
|
|
|
|
[Library subroutine to read constants; runs at load time and is then overwritten.
|
|
R5, for decimal fractions, seems to be unavailable (lost?), so the values are
|
|
here read in as 35-bit integers (i.e. times 2^34) by R2.
|
|
Values are: 0.001, initial value of y
|
|
(2^23)/(10^7) and 25/(2^10) for use in calculations
|
|
0.5/(10^9) for rounding to 9 d.p. (print routine P1 doesn't do this)]
|
|
GKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@E13Z
|
|
T#C
|
|
17179869F14411518808F419430400F9#
|
|
TZ
|
|
|
|
[Library subroutine M3; prints header at load time and is then overwritten.]
|
|
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
|
|
*SCALED!FOR!EDSAC@&!!TIME!!!!!!!!!Y!VIA!RK!!!!!Y!DIRECT@&
|
|
....PK [end text with some blank tape]
|
|
|
|
[Runge-Kutta: auxiliary subroutine to calculate (2^m)*h*(dy1/dt) and (2^m)*h*(dy2/dt)
|
|
from y1, y2, where y1 is the function y in Rosetta Code (but scaled) and y2 = t.
|
|
For the Rosetta code example we're using m = 5, h = 2^(-7)]
|
|
E25K TX GK
|
|
A3F T20@ [set up return as usual]
|
|
H2#G V2#G TD [acc := t^2, temp store in 0D]
|
|
H#G VD LD YF TD [y1 times t^2, shift left, round, temp store in 0D]
|
|
H2#C VD YF T4D [times (2^23)/(10^7), round, to 4D for square root]
|
|
[14] A14@ GV A4D T4#G [call square root, result in 4D, copy to (2^m)hy']
|
|
A21@ T6#G [1/4, i.e. (2^m)h with m and h as above, to (2^m)ht']
|
|
[20] ZF [overwritten by jump back to caller]
|
|
[21] RF [constant 1/4]
|
|
|
|
[Main routine, with two subroutines in the same address block as the main routine.]
|
|
E25K TA GK
|
|
[0] #F [figures shift on teleprinter]
|
|
[1] MF [decimal point (in figures mode)]
|
|
[2] !F @F &F [space, carriage return, line feed,]
|
|
[5] K4096F [null char]
|
|
[6] P100F [constant: nr of Runge-Kutta steps (in address field)]
|
|
[7] PF [negative count of Runge-Kutta steps]
|
|
[8] P10F [constant: number of steps between printed values]
|
|
[9] PF [negative count of steps between printed values]
|
|
[Enter with acc = 0]
|
|
[10] O@ [set teleprinter to figures]
|
|
S6@ T7@ [init negative count of R-K steps]
|
|
S8@ T9@ [init negative count of print steps]
|
|
[Before using library subroutine G1, clear its working registers (WWG page 33)]
|
|
T8#G T10#G
|
|
[Set up initial values of y1 and y2 (where y2 = t)]
|
|
A#C T#G [load 0.001 from constants section, store in y1]
|
|
T2#G [y2 = t = 0]
|
|
[20] A20@ G40@ [call subroutine to print initial values]
|
|
[Loop round Runge-Kutta steps]
|
|
[22] TF A23@ G12G [clear accumulator, call G1 for Runge-Kutta step]
|
|
A9@ A2F U9@ [update negative print count]
|
|
G33@ [skip printing if not reached 0]
|
|
S8@ T9@ [reset negative print count]
|
|
A31@ G40@ [call subroutine to print values]
|
|
[33] TF [clear accumulator]
|
|
A7@ A2F U7@ [increment negative count of Runge-Kutta steps]
|
|
G22@ [loop till count = 0]
|
|
O5@ ZF [flush teleprinter buffer; stop]
|
|
|
|
[Subroutine to print y1 as calculated (1) by Runge-Kutta (2) direct from formula]
|
|
[40] A3F T71@ [set up return as usual]
|
|
A2#G TD [latest t (= y2) from Runge-Kutta, to 0D for printing]
|
|
[44] A44@ G72@ [call subroutine to print t]
|
|
O2@ O2@ [followed by 2 spaces]
|
|
A#G TD [latest y1 from Runge-Kutta, to 0D for printing]
|
|
[50] A50@ G72@ [call subroutine to print y1]
|
|
O2@ O2@ [followed by 2 spaces]
|
|
A 4#C [load constant 25/(2^10)]
|
|
H2#G V2#G TD [add t^2, temp store result in 0D]
|
|
HD VD LD YF TD [square, shift 1 left, round, result to 0D]
|
|
H2#C VD YF TD [times (2^23)/(10^7), round, to 0D for printing]
|
|
[67] A67@ G72@ [call subroutine to print y]
|
|
O3@ O4@ [print CR, LF]
|
|
[71] ZF [overwritten by jump back to caller]
|
|
|
|
[Second-level subroutine to print number in 0D to 9 decimal places]
|
|
[72] A3F T82@ [set up return as usual]
|
|
AD A6#C TD [load number, add decimal rounding, to 0D for printing]
|
|
O81@ O1@ [print '0.' since P1 doesn't do so]
|
|
A79@ GB [call library subroutine P1 for printing]
|
|
[81] P9F [parameter for P1, 9 decimals]
|
|
[82] ZF [overwritten by jump back to caller]
|
|
|
|
[Library subroutine G1 for Runge-Kutta process. 66 locations, even address.]
|
|
E25K T12G
|
|
GKT4#ZH682DT6#ZPNT12#Z!1405DT14#ZTHT16#ZT2HTZA3FT61@A31@G63@&FT6ZPN
|
|
T8ZMMO&H4@A20@E23@T14ZAHT16ZA2HT18ZH12#@S12#@T12#@E28@H4#@T4DUFS38@
|
|
A25@T38@S6#@A16#@U46#@A8@U37@A9@U55@A24@T39@ZFR1057#@ZFYFU6DV6DRLYF
|
|
UDZFZFADLDADLLS6DN4DYFZFA46#@S14#@G29@A65@S11@ZFA35@U65@GXZF
|
|
|
|
[Replacement for library routine S2 (square root). 38 locations, even address.
|
|
Advantages: More accurate for small values of the argument.
|
|
Calculates sqrt(0) without going into an infinite loop.
|
|
Disadvantages: Longer and slower than S2 (calculates one bit at a time).]
|
|
E25K TV
|
|
GKA3FT31@A4DG32@A33@T36#@T4DA33@RDU34#@RDS4DS33@A36#@G22@T36#@A4DS34#@
|
|
T4DA36#@A33@G25@TFA36#@S33@A36#@T36#@A34#@RDYFG9@ZFZFK4096FPFPFPFPF
|
|
|
|
[Library subroutine P1 - print a single positive number. 21 locations.
|
|
Prints number in 0D to n places of decimals, where
|
|
n is specified by 'P n F' pseudo-order after subroutine call.]
|
|
E25K TB
|
|
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F
|
|
|
|
[Define entry point in main routine]
|
|
E25K TA GK
|
|
E10Z PF [enter at relative address 10 with accumulator = 0]
|