Numerical solutions to transcendental equation with nonanalytic integral

I want to find solutions y to the equation

,

for every value in C = np.linspace(0, 1000, 100). I wrote the following code but it does not work.

import numpy as np
import scipy.integrate as integrate
from scipy.optimize import root
import matplotlib.pyplot as plt

n_0 = 1
k = 5
A = np.deg2rad(30)

def n(y):
    return k*y + n_0

def f(y):
    f = lambda i: (((n(i)) / (n(0) * np.cos(A)))**2 - 1)**(-0.5)
    result, error = integrate.quad(f, 0, y)
    return result

In the above section, I have defined my function f(y). Here is the rest of my code.

C = np.linspace(0, 1000, 1000)
dy = np.array([])

for arg in C: 
    def integral(y):
        return f(y) - arg
    dy = np.append(dy, root(integral, 1).x)

plt.plot(C, dy)
plt.show()

I get the message

d:\FMF\MAF2\simulacija.py:14: RuntimeWarning: invalid value encountered in double_scalars
  f = lambda i: (((n(i)) / (n(0) * np.cos(A)))**2 - 1)**(-0.5)
d:\FMF\MAF2\simulacija.py:15: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  result, error = integrate.quad(f, 0, y)
d:\FMF\MAF2\simulacija.py:15: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  result, error = integrate.quad(f, 0, y)

Trying your code with a regular list append, I get:

In [156]: C = np.linspace(0, 1000, 1000)
...: dy = []
...: 
...: for arg in C: 
...:     def integral(y):
...:         print(y)
...:         return f(y) - arg
...:     dy.append(root(integral, 1).x)
[1]
[1.]
[1.]
[1.00000001]
[-1.46296662]
Traceback (most recent call last):
Cell In[156], line 8
dy.append(root(integral, 1).x)
File ~\miniconda3\lib\site-packages\scipy\optimize\_root.py:235 in root
sol = _root_hybr(fun, x0, args=args, jac=jac, **options)
File ~\miniconda3\lib\site-packages\scipy\optimize\_minpack_py.py:240 in _root_hybr
retval = _minpack._hybrd(func, x0, args, 1, xtol, maxfev,
Cell In[156], line 7 in integral
return f(y) - arg
Cell In[147], line 3 in f
result, error = integrate.quad(f, 0, y)
File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:463 in quad
retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:575 in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
TypeError: must be real number, not complex

So trying a specific y

In [157]: f(-1.4)
Traceback (most recent call last):
Cell In[157], line 1
f(-1.4)
Cell In[147], line 3 in f
result, error = integrate.quad(f, 0, y)
File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:463 in quad
retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:575 in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
TypeError: must be real number, not complex

quad objects if its func returns a complex.

Let's change the f def, to avoid function names over lap

In [181]: def f(y):
...:     fnc = lambda i: (((n(i)) / (n(0) * np.cos(A)))**2 - 1)**(-0.5)
...:     # print(fnc(0), fnc(y))
...:     result, error = integrate.quad(fnc, 0, y)
...:     return result
...:     
In [182]: C = np.linspace(0, 1000, 1000)
...: dy = []
...: 
...: for arg in C: 
...:     def integral(y):
...:         #print(y)
...:         return f(y) - arg
...:     dy.append(root(integral, 1).x)

Now I get some integration run time errors, and very long runs.

I haven't solved the problem, but hopefully this gives some ideas on how to debug code like this.