-1

I have a function in sympy which is quite ugly:

In [79]: print(expected_c)
Out[79]: 2**(n - 2)*n*(n - 1)*binomial(m, n)*factorial(m - 3/2)*factorial(m - n)/(binomial(2*m, n)*factorial(m - n/2)*factorial(m - n/2 - 1/2))

I want to solve it for some $m$ (say 10000) to find $n$ such that my formula values, say 4.

It is easy to do it 'by hand' by graphing it, and to see that it reaches 4 for about n=400. Then I try this value:

In [64]: expected_c.subs([(m,10000), (n,400)]).evalf()
Out[64]: 3.9901995099755

However the numerical solver fail to find this value:

In [82]: sympy.nsolve(expected_c.subs(m,10000.0)-4.0,n,400)
Out[82]: mpf('nan')

Is there a way to make it work or do I have to code a numerical solver myself ? Even a stupid one could do the job as I do not need a high precision but I was thinking Sympy should be able to do that.

Alex Riley
  • 169,130
  • 45
  • 262
  • 238
Cédric Van Rompay
  • 2,489
  • 1
  • 18
  • 24

1 Answers1

1

The docstring of nsolve warns that the numerator of your expression is used. In your case this is catastrophic and it is advisable to just use mpmath.findroot directly:

>>> m = 10**4; c = 4
>>> f=lambda n,m,c:( 2**(n - 2)*n*(n - 1)*binomial(m, n)*factorial(m - S(3)/2)*
... factorial(m - n)/(binomial(2*m, n)*factorial(m - n/2)*factorial(m - n/2
... - S(1)/2))-c)
>>> mpmath.findroot(lambda x:f(x, m, c),(300,450),solver='bisect')
mpf('400.49031238268759')
smichr
  • 16,948
  • 2
  • 27
  • 34
  • Thanks, that's great. I had tried to use findroot() before but it seems I was not using it the proper way because it was running forever so I had to stop the kernel. I was not specifying "solver='bisect'", do you think it could be the reason ? – Cédric Van Rompay Jun 03 '15 at 12:12
  • answered that myself: yes, without solver='bisect' it takes forever. – Cédric Van Rompay Jun 03 '15 at 12:13