akhhhh's profile - activity

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.

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?