Solving a complex implicit equation on python

Question

I have this function:

atanh((y/10^3)-3)-((10^(3)*y^(3)-11*10^(6)*y**2 + 3.93*10^(10)*y - 4.47*10^13)/((y - 2 * 10^3)*(y -4 * 10^3)^3))=3.49875*10^(-4)*t-atanh(3) + 0.3489583 

(As you can see in the image).

and I need it to return a y value (as a float), for a given time (t) value. I have to try the times t=0.1; t=2; t= 0,2.

I tried the fsolve function but I get an error message. This is what I did:

from cmath import atanh
from scipy.optimize import fsolve

t=1
def f(y,t):
    return atanh((y/10**3)-3)-((10**(3)*y**(3)-11*10**(6)*y**2 + 3.93*10**(10)*y - 4.47*10**13)/((y - 2 * 10**3)*(y -4 * 10**3)**3))-3.49875*10**(-4)*t+atanh(3) - 0.3489583

fsolve(f(y,t),1.9)

For which I got this error:

TypeError                                 Traceback (most recent call last)
<ipython-input-12-a340e1c537e4> in <module>
      6     return atanh((y/10**3)-3)-((10**(3)*y**(3)-11*10**(6)*y**2 + 3.93*10**(10)*y - 4.47*10**13)/((y - 2 * 10**3)*(y -4 * 10**3)**3))-3.49875*10**(-4)*t+atanh(3) - 0.3489583
      7 
----> 8 fsolve(f(y,t),1.9)

<ipython-input-12-a340e1c537e4> in f(y, t)
      4 t=1
      5 def f(y,t):
----> 6     return cmath.atanh((y/10**3)-3)-((10**(3)*y**(3)-11*10**(6)*y**2 + 3.93*10**(10)*y - 4.47*10**13)/((y - 2 * 10**3)*(y -4 * 10**3)**3))-3.49875*10**(-4)*t+cmath.atanh(3) - 0.3489583
      7 
      8 fsolve(f(y,t),1.9)

~\Anaconda3\lib\site-packages\sympy\core\expr.py in __complex__(self)
    283         result = self.evalf()
    284         re, im = result.as_real_imag()
--> 285         return complex(float(re), float(im))
    286 
    287     def __ge__(self, other):

~\Anaconda3\lib\site-packages\sympy\core\expr.py in __float__(self)
    278         if result.is_number and result.as_real_imag()[1]:
    279             raise TypeError("can't convert complex to float")
--> 280         raise TypeError("can't convert expression to float")
    281 
    282     def __complex__(self):

TypeError: can't convert expression to float

What I was hoping to get on the output was y as a real number . I did search for other similar questions on this site but I still couldn't solve this problem .

Answer 1

The sibling answer shows the way to pass the extra parameter t to f. But, then you run into another problem. The equation seems to be complex, while fsolve only works for real functions.

One way to tackle these, is mpmath, Python's multiprecision library. mpmath has a function findroot that also works for complex numbers. Note that due to the multiprecision, it can be much slower than other libraries. I don't directly see a simple way to pass the t parameter, so I made use of a lambda function:

from mpmath import findroot, atanh

def f(y,t):
    return atanh((y/10**3)-3)-((10**(3)*y**(3)-11*10**(6)*y**2 + 3.93*10**(10)*y - 4.47*10**13)/((y - 2 * 10**3)*(y -4 * 10**3)**3))-3.49875*10**(-4)*t+atanh(3) - 0.3489583

y0 = 1.9
for t in (0.1, 2, 0.2):
    ans = findroot(lambda y: f(y,t), y0)
    print(f'f({ans}, {t}) = {y0}')

output:

f((-52.8736406772712 + 1.4361714816895e-17j), 0.1) = 1.9
f((89.2356161023805 + 1.85315086887834e-19j), 2) = 1.9
f((-44.2974817249413 + 5.70332910817907e-18j), 0.2) = 1.9

I also tried to visualize the function. It seems that the real part is almost linearly dependent on t, while the imaginary part is extremely small. For some values of t, findroot doesn't find a solution within its default tolerance. I just skipped these t values. You might want to experiment with tolerances and available solvers for findroot.

Here are the code and the plot:

import numpy as np
import matplotlib.pyplot as plt

N = 10000
ts = np.linspace(0, 3, N)
real = np.zeros(N)
imag = np.zeros(N)
for i, t in enumerate(ts):
    try:
        ans = findroot(lambda y: f(y, t), y0)
    except:
        print("no solution at", t)
        pass # just use the previous value of ans
    real[i], imag[i] = ans.real, ans.imag

#plt.plot(real, imag, 'b', lw=1)
scat = plt.scatter(real, imag, c=ts, s=5)
plt.ylim((min(imag), max(imag)))
plt.xlabel('real axis')
plt.ylabel('imaginary axis')
plt.title('y values for f(y,t)=1.9, t=[0, 3]', fontsize=13)
plt.colorbar(scat, label='t values')
plt.show()