Can't find a way to solve this nonlinear equation - Printable Version +- Python Forum (https://python-forum.io) +-- Forum: Python Coding (https://python-forum.io/forum-7.html) +--- Forum: Data Science (https://python-forum.io/forum-44.html) +--- Thread: Can't find a way to solve this nonlinear equation (/thread-20512.html) |
Can't find a way to solve this nonlinear equation - Alex009988 - Aug-15-2019 I want to solve this nonlinear equation. I've tried to use sympy and scipy methods. But I thing something wrong with syntax. So how to write this code correctly? import sympy import numpy as np import math from scipy.optimize import fsolve def PoissonEqs(p,m,alpha): if m==0: return math.exp(-p)-alpha else: return math.exp(m*math.log(p)-p-math.log(math.gamma(m+1)))+PoissonEqs(p,m-1,alpha) #p=fsolve(PoissonEqs(p,0,0.05,1),1) p=sympy.Symbol('p') p1=sympy.nsolve(math.exp(-p)-0.05,p,1) print(p1) RE: Can't find a way to solve this nonlinear equation - Alex009988 - Aug-15-2019 Well seems fsolve is not a good option for this equation as I need to set a p0 guess with a certain precision. More over I have the same problem as with vba in excel when m, that is factorial, more than 170. it dosent work. It there some data types kinda long long in c++ for python? Can it be solved with nsolve? import sympy import numpy as np import math import scipy.optimize def PoissonEqs(p,m,alpha): if m==0: return math.exp(-p)-alpha else: return math.exp(m*math.log(p)-p-math.log(math.gamma(m+1)))+PoissonEqs(p,m-1,alpha) def func2(p): m=1; alpha=0.05 return PoissonEqs(p,m,alpha) #sol = scipy.optimize.bisect(func2, 0.0000000001, 1000) #sol2=scipy.optimize.fsolve(func2,10) from sympy.abc import p sol2=sympy.nsolve(math.exp(-p)-0.05,1) print(sol2) RE: Can't find a way to solve this nonlinear equation - scidam - Aug-16-2019 Not sure about how to interpret this, but at least it is better to drop using recursion in such cases; It might be quite slow, and even cause sys/recursion/limit exception. Take a look at slightly modified code, it should help you to get what you want. It seems that bisect method works as expected. import sympy import numpy as np import math from scipy.optimize import fsolve,bisect from functools import partial @np.vectorize def PoissonEqs(p, m, alpha): result = list() for j in range(m + 1): if j == 0: result.append(math.exp(-p)-alpha) continue result.append(math.exp(m*math.log(p)-p-math.log(math.gamma(m+1)))+result[-1]) return result.pop() def func2(p): m=1; alpha=0.05 return PoissonEqs(p,m,alpha) from pylab import * plot(np.linspace(0.1,10,100), func2(np.linspace(0.1,10,100)), 'r') show() bisect(func2, 0.0000000001, 10) |