On Monday, Dec. 30, 2013, I wrote:
> I find an obscure error when calling rk() inside a
> homemade function to integrate dr/dt > 0 part of
> orbit for rutherford scattering (natural units).
I have replaced the previous homemade function by
the simpler 'call_rk' function here, and focused on
why the particular value b = 0.5 gives an error when
used with call_rk, but not when rk is called interactively.
Using b = 0.49 and b = 0.51 gives no problem.
using windows xp binary, xmaxima,
maxima version 5.28.0
-------------------------------------------
Maxima 5.28.0-2 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) fpprintprec:8$
(%i2) call_rk (e, b, tmax, dt) :=
block([w, phi_inf, theta0, rmin, th, r, numer], numer:true,
w : 1/(e*b),
phi_inf : acos(w/sqrt(4 + w^2)),
theta0 : %pi - phi_inf,
rmin : b*(w + sqrt(4 + w^2))/2,
print (" theta0 = ", theta0, " rmin = ", rmin),
rk ( [-b/r^2, sqrt(1 - 1/(e*r) - b^2/r^2)], [th, r],
[theta0, rmin], [t, 0, tmax, dt]))$
------------------------------------
/* case b = 1, away from problem area */
(%i3) call_rk (1, 1, 1, 0.25);
theta0 = 2.0344439 rmin = 1.618034
(%o3) [[0,2.0344439,1.618034],[0.25,1.9389619,1.6190319],
[0.5,1.8444544,1.6384567],[0.75,1.7536003,1.6837792],
[1.0,1.6687194,1.7531982]]
---------------------------------------
/* case b = 0.49: no problem */
(%i4) call_rk (1, 0.49, 1, 0.25);
theta0 = 2.3662952 rmin = 1.2000714
(%o4) [[0,2.3662952,1.2000714],[0.25,2.2812495,1.2012438],
[0.5,2.1975403,1.2234616],[0.75,2.1186944,1.2747836],
[1.0,2.0474184,1.3525564]]
-------------------------------------
/* case b = 0.51: no problem */
(%i5) call_rk (1, 0.51, 1, 0.25);
theta0 = 2.3462938 rmin = 1.2142129
(%o5) [[0,2.3462938,1.2142129],[0.25,2.2598129,1.2142129],
[0.5,2.1733319,1.2142129],[0.75,2.0868509,1.2142129],
[1.0,2.00037,1.2142129]]
---------------------------------
/* the above cases bracket the bad case b = 0.5 */
(%i6) call_rk (1, 0.5, 1, 0.25);
theta0 = 2.3561945 rmin = 1.2071068
Expecting a number when the initial state is
replaced in the equations, but instead found:
1.05367121E-8*(-1)^0.5
#0: rk(odes = [-0.5/r^2, (-1/r - 0.25/r^2 + 1)^0.5],
state = [th, r], initial = [2.3561945, 1.2071068],
domain = [t, 0, 1, 0.25])(dynamics.mac line 392)
#1: call_rk (e = 1, b = 0.5, tmax = 1, dt = 0.25)
-- an error. To debug this try: debugmode(true);
-----------------------------------
/* go back to default fpprintprec here and repeat: */
(%i7) fpprintprec:0$
(%i8) call_rk (1, 0.5, 1, 0.25);
theta0 = 2.356194490192345 rmin = 1.207106781186548
Expecting a number when the initial state is replaced
in the equations, but instead found:
1.0536712127723509E-8*(-1)^0.5
#0: rk(odes = [-0.5/r^2, (-1/r - 0.25/r^2 + 1)^0.5],
state = [th, r],
initial = [2.356194490192345, 1.207106781186548],
domain = [t, 0, 1, 0.25])(dynamics.mac line 392)
#1: call_rk (e = 1, b = 0.5, tmax = 1, dt = 0.25)
-- an error. To debug this try: debugmode(true);
---------------------------------------------
/* calling rk interactively with same parameters gives
no problem: */
(%i9) rk ([-0.5/r^2, (-1/r - 0.25/r^2 + 1)^0.5], [th, r],
[2.356194490192345, 1.207106781186548], [t, 0, 1, 0.25]);
(%o9) [[0,2.356194490192345,1.207106781186548],
[0.25,2.270422470869927,1.20831276108977],
[0.5,2.186002526825104,1.23063846483488],
[0.75,2.106470276074839,1.282015671878286],
[1.0,2.034536801025241,1.359799098574356]]
=======================
Ted Woollett
-------------------------------
(%i13) display2d:true$
(%i14) build_info();
(%o14)
Maxima version: "5.28.0-2"
Maxima build date: "2012-08-27 23:16:48"
Host type: "i686-pc-mingw32"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.8"