ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 05 Feb 2014 08:16:57 -0600Looong computation by desolve_rk4http://ask.sagemath.org/question/10996/looong-computation-by-desolve_rk4/I'm having an issue with very long computation times for numerical solutions to DEs using desolve_rk4. It happens when computing a solution that runs off to +/- infinity. For example, the following command will do it:
t,y = var('t,y')
f(t,y) = 6+y-y^2
desolve_rk4(f(t,y), y, ivar=t, ics=[0,-4], end_points=[0,10], output='plot')
If I make "end_points" run from 0 to something much smaller, like 0.1, everything is okay. But I want some way to fix this behavior so that I can change the initial conditions easily for classroom demos without having to constantly change the endpoints (in an interact, for example). I also don't quite understand why it should be taking so long simply to run to infinity -- this behavior doesn't match my understanding of the RK4 algorithm. (I can understand that we might hit overflow quickly, but I would expect it should simply stop computation at that point.)
Any suggestions?Wed, 05 Feb 2014 07:47:15 -0600http://ask.sagemath.org/question/10996/looong-computation-by-desolve_rk4/Answer by kcrisman for <p>I'm having an issue with very long computation times for numerical solutions to DEs using desolve_rk4. It happens when computing a solution that runs off to +/- infinity. For example, the following command will do it:</p>
<pre><code>t,y = var('t,y')
f(t,y) = 6+y-y^2
desolve_rk4(f(t,y), y, ivar=t, ics=[0,-4], end_points=[0,10], output='plot')
</code></pre>
<p>If I make "end_points" run from 0 to something much smaller, like 0.1, everything is okay. But I want some way to fix this behavior so that I can change the initial conditions easily for classroom demos without having to constantly change the endpoints (in an interact, for example). I also don't quite understand why it should be taking so long simply to run to infinity -- this behavior doesn't match my understanding of the RK4 algorithm. (I can understand that we might hit overflow quickly, but I would expect it should simply stop computation at that point.) </p>
<p>Any suggestions?</p>
http://ask.sagemath.org/question/10996/looong-computation-by-desolve_rk4/?answer=16010#post-id-16010I actually get
TypeError: unable to make sense of Maxima expression ...
when I try this. In fact,
sage: desolve_rk4(f(t,y), y, ivar=t, ics=[0,-4], end_points=[0,.4])
[[0, -4],
[0.1, -6.44026308551],
[0.2, -18.4360863753],
[0.3, -9638.98851821],
[0.4, -2.29907718472e+44]]
works, but
sage: desolve_rk4(f(t,y), y, ivar=t, ics=[0,-4], end_points=[0,.5])
TypeError: unable to make sense of Maxima expression '[[0,-4],[0.1,-6.440263085506939],[0.2,-18.436086375306434],[.30000000000000004,-9638.988518207414],[0.4,-229907718471530700000000000000000000000000000.],[0.5,-.000000000000000.000000000000000]]' in Sage
and indeed `-.000000000000000.000000000000000` makes no sense inside of Sage. I'm not sure what it's supposed to mean in Maxima, either!
Here is the problem, in Maxima.
(%i1) rk(6+y-y^2,y,-4,[t,0,.5,.1]);
(%o1) [[0, - 4], [0.1, - 6.440263085506939], [0.2, - 18.436086375306434],
[.30000000000000004, - 9638.988518207414], [0.4, - 2.2990771847153071e+44],
[0.5, - .000000000000000e+4932]]
Apparently Sage doesn't know how to translate these. And what should zero times ten to the 4932 mean?
Anyway, I've opened [Trac 15789](http://trac.sagemath.org/ticket/15789) for this.Wed, 05 Feb 2014 08:16:57 -0600http://ask.sagemath.org/question/10996/looong-computation-by-desolve_rk4/?answer=16010#post-id-16010