0

I broke my problem down as follows. I am not able to solve the following equation with Python 3.9 in a meaningful way, instead it always stops with the initial_guess for small lambda_ < 1. Is there an alternative algorithm that can handle the error function better? Or can I force fsolve to search until a solution is found?

import numpy as np
from scipy.special import erfcinv, erfc
from scipy.optimize import root, fsolve

def Q(x):
    return 0.5*erfc(x/np.sqrt(2))

def Qinvers(x):
    return np.sqrt(2)*erfcinv(2*x)

def epseqn(epsilon2):
    lambda_ = 0.1
    return Q(lambda_*Qinvers(epsilon2))

eps1 = fsolve(epseqn, 1e-2)
print(eps1)

I tried root and fsolve to get a solution. Especially for the gaussian error function I do not find a solution that converges.

Myzelnetwork23
  • 3
  • 2

2 Answers2

1

root and fsolve can be used to find the roots of a function defined by f(x)=0. Since your outer function, which is basically erfc(x), has no root (it only it approaches the x-axis asymptotically from positive values) the solvers are not able to find one. Real function arguments are assumed like you did.

Mika
  • 38
  • 5
0

Before blindly starting with numerical calculations, I would recommend to think about any constraints of your function. You will find out, that your function is only defined for values between zero and one. If you assume that there is only a single root in this interval, I would recommend to use an interval search method like brentq, see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq and https://en.wikipedia.org/wiki/Brent%27s_method. However, you could instead think further and/or just plot your function, e.g. using matplotlib

import matplotlib.pyplot as plt
x = np.linspace(0, 1, 1000)
y = epseqn(x)
plt.plot(x, y)
plt.show()

There you will see that the root is at zero, which makes sense when looking at your functions, because the inverse cumulative error function is minus infinity at zero and the regular error function gives you zero at minus infinity (mathematically in the limit sense, but numerically those functions are also defined for such input values). So without any numeric calculation, you can get the root value.

Carlos Horn
  • 1,115
  • 4
  • 17