Home » Archimedes archive » Archimedes World » AW-1994-12-Disc1.adf » Disk1Dec94 » !AWDec94/Goodies/DrawBasic/!DrawBasic/Resources/Calculus/ODEs

!AWDec94/Goodies/DrawBasic/!DrawBasic/Resources/Calculus/ODEs

This website contains an archive of files for the Acorn Electron, BBC Micro, Acorn Archimedes, Commodore 16 and Commodore 64 computers, which Dominic Ford has rescued from his private collection of floppy disks and cassettes.

Some of these files were originally commercial releases in the 1980s and 1990s, but they are now widely available online. I assume that copyright over them is no longer being asserted. If you own the copyright and would like files to be removed, please contact me.

Tape/disk: Home » Archimedes archive » Archimedes World » AW-1994-12-Disc1.adf » Disk1Dec94
Filename: !AWDec94/Goodies/DrawBasic/!DrawBasic/Resources/Calculus/ODEs
Read OK: ✔
File size: 097E bytes
Load address: 0000
Exec address: 0000
Duplicates

There are 5 duplicate copies of this file in the archive:

File contents
    1REM > DrawBasic:Resources.Calculus.ODEs
    2
    3REM =======================================================
    4REM
    5REM copyright Joe Taylor June 1994
    6REM
    7REM N.B. For ERROR THROWBACK to work The first line above
    8REM      MUST have the correct pathname of this file.
    9REM      If you move the location of the file remember to
   10REM      alter this pathname.
   11REM
   12REM =======================================================
   13
   14REM ==================================================
   15REM
   16REM Solving Ordinary Differential Equations (ODEs)
   17REM using the Runge-Kutta method.
   18REM
   19REM Reference: "NUMERICAL RECIPES" CUP Press
   20REM             Flannery,Teukolsky & Vetterling
   21REM             Chapter 15.1 pp550 et. seq.
   22REM
   23REM ===================================================
   24
   25DEF FNGraphODE(function$,x,y,n%)
   26REM ------------------------------------
   27REM
   28REM Solves dy/dx=f(x,y) [=function$]
   29REM
   30REM Solution drawn passing through (x,y)
   31REM
   32REM If n%<0 then solution is drawn in
   33REM backward direction
   34REM
   35REM -----------------------------------
   36LOCAL _path
   37PROCPathBegin(_path)
   38PROC@GraphODE(function$,x,y,n%)
   39PROCPathEnd
   40=_path
   41
   42DEF PROCGraphODE(function$,x,y,n%)
   43LOCAL _path
   44PROCPathBegin(_path)
   45PROC@GraphODE(function$,x,y,n%)
   46PROCPathEnd
   47ENDPROC
   48
   49DEF PROC@GraphODE(f$,x,y,n%)
   50LOCAL a,b,u,v,h,h_max,h_next
   51PROCGetFrameCoords(a,b,u,v)
   52h_max=(b-a)/FNMaxNoofPoints
   53CASE TRUE OF
   54 WHEN n%=TRUE    : h=h_max
   55 WHEN n%=0       : h=-h_max
   56 WHEN n%=-TRUE   : h=-h_max
   57 OTHERWISE       : h=(b-a)/n%
   58ENDCASE
   59PROCMove(x,y) : dy=EVALf$
   60REPEAT
   61 X=x : Y=y : Dy=dy
   62 y+=FNRunge_Kutta(f$,x,y,h)
   63 h_next=FNRunge_KuttaStep(f$,X,Y,h)
   64 IF h>h_next THEN
   65  h=h_next
   66  y=Y+FNRunge_Kutta(f$,X,Y,h)
   67 ENDIF
   68  x+=h
   69  IF x>b THEN h=b-X : y=Y+FNRunge_Kutta(f$,X,Y,h) : x=b
   70  IF x<a THEN h=X-a : y=Y+FNRunge_Kutta(f$,X,Y,h) : x=a
   71  dy=EVALf$
   72  IF FNInFrame(X,Y) THEN PROCBezier(X+h/3,Y+h*Dy/3,x-h/3,y-h*dy/3,x,y)
   73 h=h_next
   74UNTIL x>=b OR x<=a
   75ENDPROC
   76
   77DEF FNRunge_Kutta(f$,x,y,h)
   78LOCAL k1,k2,k3,k2,Y
   79     Y=y : k1=h*EVALf$ : x+=h/2
   80y=Y+k1/2 : k2=h*EVALf$
   81y=Y+k2/2 : k3=h*EVALf$ : x+=h/2
   82y=Y+k3   : k4=h*EVALf$
   83=k1/6+k2/3+k3/3+k4/6
   84
   85DEF FNRunge_KuttaStep(f$,X,Y,h)
   86LOCAL S,delta,eps,y1,y2
   87eps=FNEpsilon
   88S=.5
   89y1=Y+FNRunge_Kutta(f$,X,Y,h/2)
   90y2=y1+FNRunge_Kutta(f$,X+h/2,y1,h/2)
   91delta=y2-y
   92IF delta=0 THEN=h
   93IF ABSdelta<=eps THEN
   94 h=S*h*ABS(eps/delta)^.2
   95 ELSE
   96 h=S*h*ABS(eps/delta)^.25
   97 IF h<0 AND x<a THEN h=x-a
   98 IF h>0 AND x>b THEN h=b-a
   99ENDIF
  100=h

)� > DrawBasic:Resources.Calculus.ODEs

=� =======================================================
�
$� copyright Joe Taylor June 1994
�
;� N.B. For ERROR THROWBACK to work The first line above
7�      MUST have the correct pathname of this file.
	;�      If you move the location of the file remember to

�      alter this pathname.
�
=� =======================================================


8� ==================================================
�
4� Solving Ordinary Differential Equations (ODEs)
#� using the Runge-Kutta method.
�
.� Reference: "NUMERICAL RECIPES" CUP Press
1�             Flannery,Teukolsky & Vetterling
-�             Chapter 15.1 pp550 et. seq.
�
9� ===================================================

!� �GraphODE(function$,x,y,n%)
*� ------------------------------------
�
&� Solves dy/dx=f(x,y) [=function$]
�
*� Solution drawn passing through (x,y)
�
 '� If n%<0 then solution is drawn in
!� backward direction
"�
#)� -----------------------------------
$� _path
%�PathBegin(_path)
& �@GraphODE(function$,x,y,n%)
'�PathEnd
(
=_path
)
*!� �GraphODE(function$,x,y,n%)
+� _path
,�PathBegin(_path)
- �@GraphODE(function$,x,y,n%)
.�PathEnd
/�
0
1� �@GraphODE(f$,x,y,n%)
2� a,b,u,v,h,h_max,h_next
3�GetFrameCoords(a,b,u,v)
4h_max=(b-a)/�MaxNoofPoints
5
Ȏ � �
6 � n%=�    : h=h_max
7 � n%=0       : h=-h_max
8 � n%=-�   : h=-h_max
9        : h=(b-a)/n%
:�
;�Move(x,y) : dy=�f$
<�
= X=x : Y=y : Dy=dy
> y+=�Runge_Kutta(f$,x,y,h)
?& h_next=�Runge_KuttaStep(f$,X,Y,h)
@ � h>h_next �
A  h=h_next
B   y=Y+�Runge_Kutta(f$,X,Y,h)
C �
D
  x+=h
E6  � x>b � h=b-X : y=Y+�Runge_Kutta(f$,X,Y,h) : x=b
F6  � x<a � h=X-a : y=Y+�Runge_Kutta(f$,X,Y,h) : x=a
G  dy=�f$
HB  � �InFrame(X,Y) � �Bezier(X+h/3,Y+h*Dy/3,x-h/3,y-h*dy/3,x,y)
I
 h=h_next
J� x>=b � x<=a
K�
L
M� �Runge_Kutta(f$,x,y,h)
N� k1,k2,k3,k2,Y
O      Y=y : k1=h*�f$ : x+=h/2
Py=Y+k1/2 : k2=h*�f$
Q y=Y+k2/2 : k3=h*�f$ : x+=h/2
Ry=Y+k3   : k4=h*�f$
S=k1/6+k2/3+k3/3+k4/6
T
U � �Runge_KuttaStep(f$,X,Y,h)
V� S,delta,eps,y1,y2
Weps=�Epsilon
XS=.5
Y!y1=Y+�Runge_Kutta(f$,X,Y,h/2)
Z'y2=y1+�Runge_Kutta(f$,X+h/2,y1,h/2)
[delta=y2-y
\� delta=0 �=h
]� �delta<=eps �
^ h=S*h*�(eps/delta)^.2
_ �
` h=S*h*�(eps/delta)^.25
a � h<0 � x<a � h=x-a
b � h>0 � x>b � h=b-a
c�
d=h
�
00000000  0d 00 01 29 f4 20 3e 20  44 72 61 77 42 61 73 69  |...). > DrawBasi|
00000010  63 3a 52 65 73 6f 75 72  63 65 73 2e 43 61 6c 63  |c:Resources.Calc|
00000020  75 6c 75 73 2e 4f 44 45  73 0d 00 02 04 0d 00 03  |ulus.ODEs.......|
00000030  3d f4 20 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 3d 3d 3d  |=. =============|
00000040  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 3d 3d 3d  |================|
*
00000060  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 0d 00 04 05 f4 0d  |==========......|
00000070  00 05 24 f4 20 63 6f 70  79 72 69 67 68 74 20 4a  |..$. copyright J|
00000080  6f 65 20 54 61 79 6c 6f  72 20 4a 75 6e 65 20 31  |oe Taylor June 1|
00000090  39 39 34 0d 00 06 05 f4  0d 00 07 3b f4 20 4e 2e  |994........;. N.|
000000a0  42 2e 20 46 6f 72 20 45  52 52 4f 52 20 54 48 52  |B. For ERROR THR|
000000b0  4f 57 42 41 43 4b 20 74  6f 20 77 6f 72 6b 20 54  |OWBACK to work T|
000000c0  68 65 20 66 69 72 73 74  20 6c 69 6e 65 20 61 62  |he first line ab|
000000d0  6f 76 65 0d 00 08 37 f4  20 20 20 20 20 20 4d 55  |ove...7.      MU|
000000e0  53 54 20 68 61 76 65 20  74 68 65 20 63 6f 72 72  |ST have the corr|
000000f0  65 63 74 20 70 61 74 68  6e 61 6d 65 20 6f 66 20  |ect pathname of |
00000100  74 68 69 73 20 66 69 6c  65 2e 0d 00 09 3b f4 20  |this file....;. |
00000110  20 20 20 20 20 49 66 20  79 6f 75 20 6d 6f 76 65  |     If you move|
00000120  20 74 68 65 20 6c 6f 63  61 74 69 6f 6e 20 6f 66  | the location of|
00000130  20 74 68 65 20 66 69 6c  65 20 72 65 6d 65 6d 62  | the file rememb|
00000140  65 72 20 74 6f 0d 00 0a  1f f4 20 20 20 20 20 20  |er to.....      |
00000150  61 6c 74 65 72 20 74 68  69 73 20 70 61 74 68 6e  |alter this pathn|
00000160  61 6d 65 2e 0d 00 0b 05  f4 0d 00 0c 3d f4 20 3d  |ame.........=. =|
00000170  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 3d 3d 3d  |================|
*
000001a0  3d 3d 3d 3d 3d 3d 0d 00  0d 04 0d 00 0e 38 f4 20  |======.......8. |
000001b0  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 3d 3d 3d  |================|
*
000001e0  3d 3d 0d 00 0f 05 f4 0d  00 10 34 f4 20 53 6f 6c  |==........4. Sol|
000001f0  76 69 6e 67 20 4f 72 64  69 6e 61 72 79 20 44 69  |ving Ordinary Di|
00000200  66 66 65 72 65 6e 74 69  61 6c 20 45 71 75 61 74  |fferential Equat|
00000210  69 6f 6e 73 20 28 4f 44  45 73 29 0d 00 11 23 f4  |ions (ODEs)...#.|
00000220  20 75 73 69 6e 67 20 74  68 65 20 52 75 6e 67 65  | using the Runge|
00000230  2d 4b 75 74 74 61 20 6d  65 74 68 6f 64 2e 0d 00  |-Kutta method...|
00000240  12 05 f4 0d 00 13 2e f4  20 52 65 66 65 72 65 6e  |........ Referen|
00000250  63 65 3a 20 22 4e 55 4d  45 52 49 43 41 4c 20 52  |ce: "NUMERICAL R|
00000260  45 43 49 50 45 53 22 20  43 55 50 20 50 72 65 73  |ECIPES" CUP Pres|
00000270  73 0d 00 14 31 f4 20 20  20 20 20 20 20 20 20 20  |s...1.          |
00000280  20 20 20 46 6c 61 6e 6e  65 72 79 2c 54 65 75 6b  |   Flannery,Teuk|
00000290  6f 6c 73 6b 79 20 26 20  56 65 74 74 65 72 6c 69  |olsky & Vetterli|
000002a0  6e 67 0d 00 15 2d f4 20  20 20 20 20 20 20 20 20  |ng...-.         |
000002b0  20 20 20 20 43 68 61 70  74 65 72 20 31 35 2e 31  |    Chapter 15.1|
000002c0  20 70 70 35 35 30 20 65  74 2e 20 73 65 71 2e 0d  | pp550 et. seq..|
000002d0  00 16 05 f4 0d 00 17 39  f4 20 3d 3d 3d 3d 3d 3d  |.......9. ======|
000002e0  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 3d 3d 3d  |================|
*
00000300  3d 3d 3d 3d 3d 3d 3d 3d  3d 3d 3d 3d 3d 0d 00 18  |=============...|
00000310  04 0d 00 19 21 dd 20 a4  47 72 61 70 68 4f 44 45  |....!. .GraphODE|
00000320  28 66 75 6e 63 74 69 6f  6e 24 2c 78 2c 79 2c 6e  |(function$,x,y,n|
00000330  25 29 0d 00 1a 2a f4 20  2d 2d 2d 2d 2d 2d 2d 2d  |%)...*. --------|
00000340  2d 2d 2d 2d 2d 2d 2d 2d  2d 2d 2d 2d 2d 2d 2d 2d  |----------------|
00000350  2d 2d 2d 2d 2d 2d 2d 2d  2d 2d 2d 2d 0d 00 1b 05  |------------....|
00000360  f4 0d 00 1c 26 f4 20 53  6f 6c 76 65 73 20 64 79  |....&. Solves dy|
00000370  2f 64 78 3d 66 28 78 2c  79 29 20 5b 3d 66 75 6e  |/dx=f(x,y) [=fun|
00000380  63 74 69 6f 6e 24 5d 0d  00 1d 05 f4 0d 00 1e 2a  |ction$]........*|
00000390  f4 20 53 6f 6c 75 74 69  6f 6e 20 64 72 61 77 6e  |. Solution drawn|
000003a0  20 70 61 73 73 69 6e 67  20 74 68 72 6f 75 67 68  | passing through|
000003b0  20 28 78 2c 79 29 0d 00  1f 05 f4 0d 00 20 27 f4  | (x,y)....... '.|
000003c0  20 49 66 20 6e 25 3c 30  20 74 68 65 6e 20 73 6f  | If n%<0 then so|
000003d0  6c 75 74 69 6f 6e 20 69  73 20 64 72 61 77 6e 20  |lution is drawn |
000003e0  69 6e 0d 00 21 18 f4 20  62 61 63 6b 77 61 72 64  |in..!.. backward|
000003f0  20 64 69 72 65 63 74 69  6f 6e 0d 00 22 05 f4 0d  | direction.."...|
00000400  00 23 29 f4 20 2d 2d 2d  2d 2d 2d 2d 2d 2d 2d 2d  |.#). -----------|
00000410  2d 2d 2d 2d 2d 2d 2d 2d  2d 2d 2d 2d 2d 2d 2d 2d  |----------------|
00000420  2d 2d 2d 2d 2d 2d 2d 2d  0d 00 24 0b ea 20 5f 70  |--------..$.. _p|
00000430  61 74 68 0d 00 25 15 f2  50 61 74 68 42 65 67 69  |ath..%..PathBegi|
00000440  6e 28 5f 70 61 74 68 29  0d 00 26 20 f2 40 47 72  |n(_path)..& .@Gr|
00000450  61 70 68 4f 44 45 28 66  75 6e 63 74 69 6f 6e 24  |aphODE(function$|
00000460  2c 78 2c 79 2c 6e 25 29  0d 00 27 0c f2 50 61 74  |,x,y,n%)..'..Pat|
00000470  68 45 6e 64 0d 00 28 0a  3d 5f 70 61 74 68 0d 00  |hEnd..(.=_path..|
00000480  29 04 0d 00 2a 21 dd 20  f2 47 72 61 70 68 4f 44  |)...*!. .GraphOD|
00000490  45 28 66 75 6e 63 74 69  6f 6e 24 2c 78 2c 79 2c  |E(function$,x,y,|
000004a0  6e 25 29 0d 00 2b 0b ea  20 5f 70 61 74 68 0d 00  |n%)..+.. _path..|
000004b0  2c 15 f2 50 61 74 68 42  65 67 69 6e 28 5f 70 61  |,..PathBegin(_pa|
000004c0  74 68 29 0d 00 2d 20 f2  40 47 72 61 70 68 4f 44  |th)..- .@GraphOD|
000004d0  45 28 66 75 6e 63 74 69  6f 6e 24 2c 78 2c 79 2c  |E(function$,x,y,|
000004e0  6e 25 29 0d 00 2e 0c f2  50 61 74 68 45 6e 64 0d  |n%).....PathEnd.|
000004f0  00 2f 05 e1 0d 00 30 04  0d 00 31 1b dd 20 f2 40  |./....0...1.. .@|
00000500  47 72 61 70 68 4f 44 45  28 66 24 2c 78 2c 79 2c  |GraphODE(f$,x,y,|
00000510  6e 25 29 0d 00 32 1c ea  20 61 2c 62 2c 75 2c 76  |n%)..2.. a,b,u,v|
00000520  2c 68 2c 68 5f 6d 61 78  2c 68 5f 6e 65 78 74 0d  |,h,h_max,h_next.|
00000530  00 33 1c f2 47 65 74 46  72 61 6d 65 43 6f 6f 72  |.3..GetFrameCoor|
00000540  64 73 28 61 2c 62 2c 75  2c 76 29 0d 00 34 1e 68  |ds(a,b,u,v)..4.h|
00000550  5f 6d 61 78 3d 28 62 2d  61 29 2f a4 4d 61 78 4e  |_max=(b-a)/.MaxN|
00000560  6f 6f 66 50 6f 69 6e 74  73 0d 00 35 0a c8 8e 20  |oofPoints..5... |
00000570  b9 20 ca 0d 00 36 18 20  c9 20 6e 25 3d b9 20 20  |. ...6. . n%=.  |
00000580  20 20 3a 20 68 3d 68 5f  6d 61 78 0d 00 37 1c 20  |  : h=h_max..7. |
00000590  c9 20 6e 25 3d 30 20 20  20 20 20 20 20 3a 20 68  |. n%=0       : h|
000005a0  3d 2d 68 5f 6d 61 78 0d  00 38 19 20 c9 20 6e 25  |=-h_max..8. . n%|
000005b0  3d 2d b9 20 20 20 3a 20  68 3d 2d 68 5f 6d 61 78  |=-.   : h=-h_max|
000005c0  0d 00 39 19 20 7f 20 20  20 20 20 20 20 3a 20 68  |..9. .       : h|
000005d0  3d 28 62 2d 61 29 2f 6e  25 0d 00 3a 05 cb 0d 00  |=(b-a)/n%..:....|
000005e0  3b 17 f2 4d 6f 76 65 28  78 2c 79 29 20 3a 20 64  |;..Move(x,y) : d|
000005f0  79 3d a0 66 24 0d 00 3c  05 f5 0d 00 3d 16 20 58  |y=.f$..<....=. X|
00000600  3d 78 20 3a 20 59 3d 79  20 3a 20 44 79 3d 64 79  |=x : Y=y : Dy=dy|
00000610  0d 00 3e 1e 20 79 2b 3d  a4 52 75 6e 67 65 5f 4b  |..>. y+=.Runge_K|
00000620  75 74 74 61 28 66 24 2c  78 2c 79 2c 68 29 0d 00  |utta(f$,x,y,h)..|
00000630  3f 26 20 68 5f 6e 65 78  74 3d a4 52 75 6e 67 65  |?& h_next=.Runge|
00000640  5f 4b 75 74 74 61 53 74  65 70 28 66 24 2c 58 2c  |_KuttaStep(f$,X,|
00000650  59 2c 68 29 0d 00 40 11  20 e7 20 68 3e 68 5f 6e  |Y,h)..@. . h>h_n|
00000660  65 78 74 20 8c 0d 00 41  0e 20 20 68 3d 68 5f 6e  |ext ...A.  h=h_n|
00000670  65 78 74 0d 00 42 20 20  20 79 3d 59 2b a4 52 75  |ext..B   y=Y+.Ru|
00000680  6e 67 65 5f 4b 75 74 74  61 28 66 24 2c 58 2c 59  |nge_Kutta(f$,X,Y|
00000690  2c 68 29 0d 00 43 06 20  cd 0d 00 44 0a 20 20 78  |,h)..C. ...D.  x|
000006a0  2b 3d 68 0d 00 45 36 20  20 e7 20 78 3e 62 20 8c  |+=h..E6  . x>b .|
000006b0  20 68 3d 62 2d 58 20 3a  20 79 3d 59 2b a4 52 75  | h=b-X : y=Y+.Ru|
000006c0  6e 67 65 5f 4b 75 74 74  61 28 66 24 2c 58 2c 59  |nge_Kutta(f$,X,Y|
000006d0  2c 68 29 20 3a 20 78 3d  62 0d 00 46 36 20 20 e7  |,h) : x=b..F6  .|
000006e0  20 78 3c 61 20 8c 20 68  3d 58 2d 61 20 3a 20 79  | x<a . h=X-a : y|
000006f0  3d 59 2b a4 52 75 6e 67  65 5f 4b 75 74 74 61 28  |=Y+.Runge_Kutta(|
00000700  66 24 2c 58 2c 59 2c 68  29 20 3a 20 78 3d 61 0d  |f$,X,Y,h) : x=a.|
00000710  00 47 0c 20 20 64 79 3d  a0 66 24 0d 00 48 42 20  |.G.  dy=.f$..HB |
00000720  20 e7 20 a4 49 6e 46 72  61 6d 65 28 58 2c 59 29  | . .InFrame(X,Y)|
00000730  20 8c 20 f2 42 65 7a 69  65 72 28 58 2b 68 2f 33  | . .Bezier(X+h/3|
00000740  2c 59 2b 68 2a 44 79 2f  33 2c 78 2d 68 2f 33 2c  |,Y+h*Dy/3,x-h/3,|
00000750  79 2d 68 2a 64 79 2f 33  2c 78 2c 79 29 0d 00 49  |y-h*dy/3,x,y)..I|
00000760  0d 20 68 3d 68 5f 6e 65  78 74 0d 00 4a 11 fd 20  |. h=h_next..J.. |
00000770  78 3e 3d 62 20 84 20 78  3c 3d 61 0d 00 4b 05 e1  |x>=b . x<=a..K..|
00000780  0d 00 4c 04 0d 00 4d 1c  dd 20 a4 52 75 6e 67 65  |..L...M.. .Runge|
00000790  5f 4b 75 74 74 61 28 66  24 2c 78 2c 79 2c 68 29  |_Kutta(f$,x,y,h)|
000007a0  0d 00 4e 13 ea 20 6b 31  2c 6b 32 2c 6b 33 2c 6b  |..N.. k1,k2,k3,k|
000007b0  32 2c 59 0d 00 4f 20 20  20 20 20 20 59 3d 79 20  |2,Y..O      Y=y |
000007c0  3a 20 6b 31 3d 68 2a a0  66 24 20 3a 20 78 2b 3d  |: k1=h*.f$ : x+=|
000007d0  68 2f 32 0d 00 50 17 79  3d 59 2b 6b 31 2f 32 20  |h/2..P.y=Y+k1/2 |
000007e0  3a 20 6b 32 3d 68 2a a0  66 24 0d 00 51 20 79 3d  |: k2=h*.f$..Q y=|
000007f0  59 2b 6b 32 2f 32 20 3a  20 6b 33 3d 68 2a a0 66  |Y+k2/2 : k3=h*.f|
00000800  24 20 3a 20 78 2b 3d 68  2f 32 0d 00 52 17 79 3d  |$ : x+=h/2..R.y=|
00000810  59 2b 6b 33 20 20 20 3a  20 6b 34 3d 68 2a a0 66  |Y+k3   : k4=h*.f|
00000820  24 0d 00 53 18 3d 6b 31  2f 36 2b 6b 32 2f 33 2b  |$..S.=k1/6+k2/3+|
00000830  6b 33 2f 33 2b 6b 34 2f  36 0d 00 54 04 0d 00 55  |k3/3+k4/6..T...U|
00000840  20 dd 20 a4 52 75 6e 67  65 5f 4b 75 74 74 61 53  | . .Runge_KuttaS|
00000850  74 65 70 28 66 24 2c 58  2c 59 2c 68 29 0d 00 56  |tep(f$,X,Y,h)..V|
00000860  17 ea 20 53 2c 64 65 6c  74 61 2c 65 70 73 2c 79  |.. S,delta,eps,y|
00000870  31 2c 79 32 0d 00 57 10  65 70 73 3d a4 45 70 73  |1,y2..W.eps=.Eps|
00000880  69 6c 6f 6e 0d 00 58 08  53 3d 2e 35 0d 00 59 21  |ilon..X.S=.5..Y!|
00000890  79 31 3d 59 2b a4 52 75  6e 67 65 5f 4b 75 74 74  |y1=Y+.Runge_Kutt|
000008a0  61 28 66 24 2c 58 2c 59  2c 68 2f 32 29 0d 00 5a  |a(f$,X,Y,h/2)..Z|
000008b0  27 79 32 3d 79 31 2b a4  52 75 6e 67 65 5f 4b 75  |'y2=y1+.Runge_Ku|
000008c0  74 74 61 28 66 24 2c 58  2b 68 2f 32 2c 79 31 2c  |tta(f$,X+h/2,y1,|
000008d0  68 2f 32 29 0d 00 5b 0e  64 65 6c 74 61 3d 79 32  |h/2)..[.delta=y2|
000008e0  2d 79 0d 00 5c 11 e7 20  64 65 6c 74 61 3d 30 20  |-y..\.. delta=0 |
000008f0  8c 3d 68 0d 00 5d 13 e7  20 94 64 65 6c 74 61 3c  |.=h..].. .delta<|
00000900  3d 65 70 73 20 8c 0d 00  5e 1a 20 68 3d 53 2a 68  |=eps ...^. h=S*h|
00000910  2a 94 28 65 70 73 2f 64  65 6c 74 61 29 5e 2e 32  |*.(eps/delta)^.2|
00000920  0d 00 5f 06 20 cc 0d 00  60 1b 20 68 3d 53 2a 68  |.._. ...`. h=S*h|
00000930  2a 94 28 65 70 73 2f 64  65 6c 74 61 29 5e 2e 32  |*.(eps/delta)^.2|
00000940  35 0d 00 61 18 20 e7 20  68 3c 30 20 80 20 78 3c  |5..a. . h<0 . x<|
00000950  61 20 8c 20 68 3d 78 2d  61 0d 00 62 18 20 e7 20  |a . h=x-a..b. . |
00000960  68 3e 30 20 80 20 78 3e  62 20 8c 20 68 3d 62 2d  |h>0 . x>b . h=b-|
00000970  61 0d 00 63 05 cd 0d 00  64 06 3d 68 0d ff        |a..c....d.=h..|
0000097e

© Dominic Ford 2021–2025.

Our privacy policy is here.
Last updated: 01 May 2025, 21:49 UTC
Website designed by .