Ask Your Question

Looong computation by desolve_rk4

asked 2014-02-05 07:47:15 -0500

Andy gravatar image

updated 2014-02-06 04:03:26 -0500

tmonteil gravatar image

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?

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted

answered 2014-02-05 08:16:57 -0500

kcrisman gravatar image

I 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 for this.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools



Asked: 2014-02-05 07:47:15 -0500

Seen: 100 times

Last updated: Feb 05 '14