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()

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)

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
尝试特定的 `y` 值：
```python
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` 对象如果其 `func` 返回复数会出错。
让我们修改函数 `f` 的定义，以避免函数名称重叠：
```python
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)