Ask Your Question

akhhhh's profile - activity

2021-01-15 22:54:19 +0100 received badge  Student (source)
2021-01-15 19:54:18 +0100 commented answer find_root fails in for loop with type coercion error

Thanks so much Emmanuel! And thank you for the suggestion about using try and except constructs for the instances with no root--as implemented the model is actually incomplete, so hopefully when I finish setting it up it will have roots everywhere as expected, but in the mean time the way I have it set up should allow me to get a better understanding of its behavior before building it out fully.

2021-01-15 19:48:29 +0100 received badge  Scholar (source)
2021-01-15 18:13:02 +0100 asked a question find_root fails in for loop with type coercion error

Hi there, I'm trying to reproduce a model in the literature. Here is the code that causes the error (description+stack trace follows):

electron_chg = 1.6*10**(-19) # Coulombs
eps0 = 8.85418782e-12 # m^-3 kg^-1 s^4 A^2 = F/m
angs_per_meter = 1*10**10
T1 = 108
T11 = 120
T111 = 120
a10 = 3.34e5
a110 = 4.69e6
a1110 = -5.52e7
a12 = 3.23e8
a112 = 4.47e9
a123 = 4.91e9
m110 = -202e6
m1110 = 276e7

var('Px,Py,Pz, T')
a1(T) = a10*(T - T1)
a11(T) = a110*(T - T11)+m110
a111(T) = a1110*(T - T111)+m1110
Glan(Px,Py,Pz, T) = (a1(T))*(Px**2+Py**2+Pz**2) + (a11(T))*(Px**4+Py**4+Pz**4) + (a111(T))*(Px**6+Py**6+Pz**6) + a12*((Px**2)*(Py**2)+(Pz**2)*(Py**2)+(Px**2)*(Pz**2)) + a112*((Px**2)*(Py**4+Pz**4)+(Py**2)*(Px**4+Pz**4)+(Pz**2)*(Px**4+Py**4)) + a123*(Px**2)*(Py**2)*(Pz**2)

Ex(Px,Py,Pz, T) = derivative(Glan,Px)
Ey(Px,Py,Pz, T) = derivative(Glan,Py)
Ez(Px,Py,Pz, T) = derivative(Glan,Pz)

# this works
var('Ptet')
find_root(Ez(0,0,Ptet,32),0.0001,1.0)

# this doesn't, even for i = 0, which is the case that works above
var('Ptet')
Ptetx = zero_vector(101)
Ptety = zero_vector(101)
Ptetz = zero_vector(101)
for i in range(101):
    print(i)
    Ptetx[i] = find_root(Ex(0,0,Ptet,32+i),0.0001,1.0)
    Ptety[i] = find_root(Ey(0,0,Ptet,32+i),0.0001,1.0)
    Ptetz[i] = find_root(Ez(0,0,Ptet,32+i),0.0001,1.0)

Ptetx

Here is the full trace of the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-9-3d8c5561fbc3> in <module>
      5 for i in range(Integer(101)):
      6     print(i)
----> 7     Ptetx[i] = find_root(Ex(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
      8     Ptety[i] = find_root(Ey(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
      9     Ptetz[i] = find_root(Ez(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/modules/free_module_element.pyx in sage.modules.free_module_element.FreeModuleElement.__setitem__ (build/cythonized/sage/modules/free_module_element.c:13856)()
   1813             if n < 0 or n >= d:
   1814                 raise IndexError("vector index out of range")
-> 1815             self.set_unsafe(n, R(value))
   1816 
   1817     cdef int set_unsafe(self, Py_ssize_t i, value) except -1:

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9337)()
    898         if mor is not None:
    899             if no_extra_args:
--> 900                 return mor._call_(x)
    901             else:
    902                 return mor._call_with_args(x, args, kwds)

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6043)()
    285             raise TypeError("Cannot coerce {} to {}".format(x, C))
    286         cdef Map m
--> 287         cdef Element e = method(C)
    288         if e is None:
    289             raise RuntimeError("BUG in coercion model: {} method of {} returned None".format(self.method_name, type(x)))

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/rings/real_mpfr.pyx in sage.rings.real_mpfr.RealNumber._integer_ (build/cythonized/sage/rings/real_mpfr.c:17319)()
   2227             return n
   2228 
-> 2229         raise TypeError("Attempt to coerce non-integral RealNumber to Integer")
   2230 
   2231     def integer_part(self):

TypeError: Attempt to coerce non-integral RealNumber to Integer

The error describes itself as a type coercion, but I can't see what the origin is, given that the 0 case works outside the loop. Could it be related to how range() works, is it different than the python implementation?