1

So I have this complicated equation which I need to solve. I think that finally x should be of order 1E22. But the problem with this code is that it crashes my entire system. Is there a fix? I tried scipy.optimize.root but it doesn't really solve anything at this order of magnitude (it gives final answer as initial guess without any iteration).

from scipy.optimize import fsolve
import math
import mpmath
import scipy
import sympy
from sympy.solvers import solve
from sympy import Symbol
from sympy import sqrt,exp

x = Symbol('x',positive=True)

cs = 507.643E-12
esi = 1.05E-10 
q = 1.6E-19
T = 300
k = 1.381E-23
ni = 1.45E16

print(solve(exp(x/((2*cs/(esi*q))**2)) - ((x/ni)**(esi*k*T)),x))

def func(N):
    return (math.exp(N/math.pow(2*cs/(esi*q),2)) - math.pow(N/ni,(esi*k*T)))

n_initial_guess = 1E21
n_solution = fsolve(func, n_initial_guess)

print ("The solution is n = %f" % n_solution)
print ("at which the value of the expression is %f" % func(n_solution))
print(scipy.optimize.root(func, 1E22,tol=1E-10))

Neither of the scipy functions work. The sympy function crashes my laptop. Would Matlab be ideal for this?

Black Jack 21
  • 315
  • 4
  • 19
  • The problem is here `print(solve(exp(x/((2*cs/(esi*q))**2)) - ((x/ni)**(esi*k*T)),x))`. this is not the optimal way to get the result. it will basically consume a huge chunk of memory making your laptop to crash. try using a different technique to solve the equation rather than using sympy solve. – Khalil Al Hooti Oct 13 '18 at 05:37

1 Answers1

3

Numeric solution with SciPy

The problem that SciPy has with this equation is the loss of significance. You are raising N to the tiny power esi*k*T which makes it very near 1; in floating point arithmetics, it becomes exactly 1. Similarly, the part coming from the exponential becomes 1. Then the two parts are subtracting, leaving 0 - equation appears to be already solved. You could have seen this happening by printing func(1E21) -- it returns 0.

The way to deal with the loss of significance is to rewrite the equation, from the original form 

exp(x/((2*cs/(esi*q))**2)) == (x/ni)**(esi*k*T)

by raising both sides to the power 1/(esi*k*T): 

exp(x*esi*q**2/(2*cs*k*T)**2)) == x/ni

So func becomes

def func(N):
    return np.exp(N*esi*q**2/(k*T*(2*cs)**2)) - (N/ni)

(It's is advisable to use NumPy functions with SciPy solvers.) That said, the solvers, for example root(func, 1E10), will report being unable to converge to a solution.

Symbolic solution with SymPy

SymPy is for solving equations analytically. It does not care for a bunch of floating point numbers. Give it a symbolic equation instead: 

x, a, b, c = symbols('x, a, b, c', positive=True)
sol = solve(exp(x/a) - (x/b)**c, x)[0]

The solution is obtained as -c*LambertW(-a/(b*c))/a. Then it can be evaluated: 

cs = 507.643E-12
esi = 1.05E-10 
q = 1.6E-19
T = 300
k = 1.381E-23
ni = 1.45E16

print(sol.evalf(subs={a: (2*cs/(esi*q))**2, b: ni, c: esi*k*T}))

Which prints -21301663061.0653 - 4649834682.69762*I confirming what one would already expect from the failure of convergence with SciPy: there are no real solutions of the equation.