Bottom Page

Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
 Can't find a way to solve this nonlinear equation
#1
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)
Quote
#2
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)
Quote
#3
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)
Quote

Top Page

Possibly Related Threads...
Thread Author Replies Views Last Post
  odeint to solve Schrodinger equation kiyoshi7 14 2,176 Nov-23-2018, 11:49 AM
Last Post: kiyoshi7

Forum Jump:


Users browsing this thread: 1 Guest(s)