search

upiterbarg / mpmath

Automatically exported from code.google.com/p/mpmath
Other
0 stars 0 forks source link

odefun stalls #188

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago 
Attempting to solve the following system
>>> def sys(x, y):
     mu = mpf("0.012277471")
     mup = 1 - mu           

     d1 = ((y[0] + mu)**2 + (y[2])**2)**(3.0/2.0)
     d2 = ((y[0] - mup)**2 + (y[2])**2)**(3.0/2.0)

     return [y[1], y[0] + 2*y[3] - mup*(y[0] + mu) / d1 - mu*(y[0] - mup) /
d2, y[3], y[2] - 2*y[1] - mup*(y[2]) / d1 - mu*(y[2]) / d2]

as follows:
>>> f = odefun(sys, 0, [0, 0, 0, mpf("-2.00158510637908252240537862224")])
>>> f(mpf("17.0652165601579625588917206249"))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mpmath/calculus/odes.py", line 273, in interpolant
    ser, xa, xb = get_series(x)
  File "mpmath/calculus/odes.py", line 262, in get_series
    ser, xb = ode_taylor(ctx, F, xb, y, tol_prec, degree)
  File "mpmath/calculus/odes.py", line 19, in ode_taylor
    fxy = derivs(x, y)
  File "<stdin>", line 8, in sys
  File "<string>", line 7, in __mul__
  File "mpmath/libmp/libmpf.py", line 796, in python_mpf_mul
    man = sman*tman

and was left for an hour on a 2.8Ghz 64-bit Linux system but did not
return; attempted at both 50 and 6 dps; memory usage hit 1GiB.

I am unsure what properties of the problem are responsible for this,
usually odefun performs quite well for me. (On a side note: the above
system can be solved at that point using 8000 or so RK4 steps to give a
reasonably accurate answer.)

Original issue reported on code.google.com by fred...@witherden.org on 16 Mar 2010 at 11:20

GoogleCodeExporter commented 9 years ago 
Thanks for the bug report!

This could be a limitation of the Taylor series algorithm, or just a problem 
with the
implementation. Perhaps it's due to the fact that (y[2])**2)**(3.0/2.0) isn't
analytic at the starting point y[2] = 0. odefun has very little 
debugging/tweaking
support, so it probably needs some modifications just to investigate this.

I don't really have time to look at this currently, so help would be 
appreciated.

Original comment by fredrik....@gmail.com on 16 Mar 2010 at 11:51