Runge-Kutta method: Difference between revisions
Content added Content deleted
(Added solution for EDSAC.) |
|||
Line 610: | Line 610: | ||
y(9.00) = 451.56245928 Error = 9.0182772312e-8 |
y(9.00) = 451.56245928 Error = 9.0182772312e-8 |
||
y(10.0) = 675.99994902 Error = 7.5419063100e-8 |
y(10.0) = 675.99994902 Error = 7.5419063100e-8 |
||
</pre> |
|||
=={{header|EDSAC order code}}== |
|||
The EDSAC subroutine library had two Runge-Kutta subroutines: G1 for 35-bit values and G2 for 17-bit values. A demo of G1 is given here. Setting up the parameters is rather complicated, but after that it's just a matter of calling G1 once for every step in the Runge-Kutta process. |
|||
Since EDSAC real numbers are restricted to -1 <= x < 1, the values in the Rosetta Code task have to be scaled down. For comparison with other languages it's convenient to divide the y values by 1000. With 100 steps, a convenient time interval is 1/128. |
|||
G1 can solve equations in several variables, say y_1, ..., y_n. The user must provide an auxiliary subroutine which calculates dy_1/dt, ..., dy_n/dt from y_1, ..., y_n. If the derivatives also depend on t (as in the Rosetta Code task) it's necessary to add a dummy y variable which is identical with t. |
|||
<lang edsac> |
|||
[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' are stored in n consecutive long locations, the last of which is bD. |
|||
G1 uses working variables stored 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. Choose a base address, 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. (EDSAC 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: number 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@&FT6ZPNT8ZMMO&H4@A20@E23@ |
|||
T14ZAHT16ZA2HT18ZH12#@S12#@T12#@E28@H4#@T4DUFS38@A25@T38@S6#@A16#@U46#@A8@U37@A9@U55@ |
|||
A24@T39@ZFR1057#@ZFYFU6DV6DRLYFUDZFZFADLDADLLS6DN4DYFZFA46#@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] |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
SCALED FOR EDSAC |
|||
TIME Y VIA RK Y DIRECT |
|||
0.000000000 0.001000000 0.001000000 |
|||
0.078125000 0.001562499 0.001562500 |
|||
0.156250000 0.003999998 0.004000000 |
|||
0.234375000 0.010562495 0.010562500 |
|||
0.312500000 0.024999992 0.025000000 |
|||
0.390625000 0.052562487 0.052562500 |
|||
0.468750000 0.099999981 0.100000000 |
|||
0.546875000 0.175562474 0.175562500 |
|||
0.625000000 0.288999965 0.289000000 |
|||
0.703125000 0.451562456 0.451562500 |
|||
0.781250000 0.675999945 0.676000000 |
|||
</pre> |
</pre> |
||