Runge-Kutta method: Difference between revisions
Added solution for EDSAC.
(Added solution for EDSAC.) |
|||
Line 610:
y(9.00) = 451.56245928 Error = 9.0182772312e-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>
| |||