Python Forum
Can't find a way to solve this nonlinear equation
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)
Reply
#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)
Reply
#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)
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  phase portrait - Nonlinear damped pendulum Giovanni_62 1 822 Sep-10-2023, 02:25 PM
Last Post: deanhystad
  Odeint to solve Mukhanov equation Messier087 4 2,508 Feb-17-2020, 05:05 PM
Last Post: Messier087
  odeint to solve Schrodinger equation kiyoshi7 14 12,105 Nov-23-2018, 11:49 AM
Last Post: kiyoshi7

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020