"Polynomial division failed" issue

[Carlos Bouthelier Madre]

Hi,

 I've got something like f(Beta,eps)=g(  integral(h(Beta,eps,tita1,tita2,...titaN)dtita1*dtita2*....*dtitaN),  Beta,  eps) where g and h are functions, and the integral of h showed as a variable of g has no primitive.

 I want to plot f, with Beta and eps in de x and y axes. Now I am just trying to plot with eps fixed and  Beta in the X axis, but I get "PolynomialDivisionFailed" when it comes to plotting, I have followed "the calls" trying to ferret out how to fix it, but I am kind of a beginner and I'm a bit overwhelmed. I would really appreciate your help. 

Thanks in advance! Have a nice day:)

This is my code:

from numpy import *

from sympy import *

from sympy import cos

from sympy.tensor.array import Array

from sympy import I, Matrix, symbols

from sympy import init_printing

from sympy import init_session

from sympy import plot

from IPython.display import display, Math, Latex

var('Tita1 Tita2 Tita3 Tita4 Tita5 Beta eps', real=True)

def delta(n):

    f='Tita'+str(n)

    return sqrt(cos(f)*cos(f)*eps*eps+1)

def EF(i):

    return Matrix([[0,0],[0,2*delta(i)]])

ID=Matrix([[1,0],[0,1]])

def ExtDim(X,n):

    if n==1:

        return X

    else:

        return kron(ExtDim(X,n-1),ID)

def n_energy(n):

    if n==1:

        return EF(n)

    else:

        Etemp=kron(ID,n_energy(n-1))

        Etemp+=ExtDim(EF(n),n)

        return Etemp

def AsignaDiag(E,n):

    for i in range(2**n):

        Evec[i]=simplify(E[i,i]-E[2**n-1,2**n-1]/2)

    

def Z_Q(n): 

    Zq=0

    

    for i in range(2**n):

        prod=1

        for j in range(2**n):

            if j!=i:

                prod=prod*(Beta*Evec[j]-Beta*Evec[i]);

        prod=simplify(prod)

        Zq=Zq+exp(-Beta*Evec[i])/prod

    Zq=Zq*((2*pi)**(2**n))

    return simplify(Zq)

    

def Z_T(n):

    Z=Z_Q(n)

    for i in range (1,n+1):

        Z=Integral(Z,('Tita'+str(i),0,2*pi))

    return simplify(Z/Beta**(n))

def Rho_Q(n):

    x=[0]*(2**n)

    Zt=Z_T(n)

    for i in range (2**n):

        for k in range (2**n):

            Prod=1

            for j in range(2**n):

                if (j!=k):

                    Prod=Prod*(Evec[j]-Evec[k])

            Prod=simplify(Prod)

            if (k!=i):

                x[i]+=exp(-Beta*Evec[k])/((Evec[i]-Evec[k])*Prod)

            else:

                Sumprodi=0

                for l in range(2**n):

                     Prodi=1

                     for j in range(2**n):

                        if ((j!=l)and(j!=i)):

                            Prodi=Prodi*(Evec[j]-Evec[i])

                     Sumprodi+=simplify(Prodi)

                Sumprodi=simplify(Sumprodi)

                x[i]+=Beta*exp(-Beta*Evec[k])/Prod-exp(-Beta*Evec[i])*Sumprodi/(Prod**2)

        x[i]+=(2*pi/Beta)**(2**n)/Zt

    return x

def Rho_T(n):

    x=[None]*(2**n)

    x=Rho_Q(n)

    for i in range (2**n):

        for j in range (1,n+1):

            x[i]=(integral(x[i],("Tita"+str(j),0,2*pi)))

        x[i]=simplify((x[i]/(Beta**n))/Z_T(n))

    return x

def S_VN(n,x): #with Rho_T =  S_quant|Trace out clas

    S=0

    for i in range (2**n):

        S+=x[i]*log(x[i])

    return S

def S_T(n):

    ro=Rho_Q(n)

    Sf=S_VN(n,ro)

    for i in range(1,n+1):

        Sf=(Integral(Sf,('Tita'+str(i),0,2*pi)))

    Sf=(Sf/(Beta**n))

    return Sf

    

init_session()

N=2

E=n_energy(N)

init_printing()

Evec=[None]*(2**N)

AsignaDiag(E,N)

S=S_T(N)

SE=S.subs(eps,0.5)

plot(SE, (Beta,0.1,1000))#This is the plot I want to make

[Aaron Meurer]

On Wed, May 10, 2017 at 2:19 PM, Carlos Bouthelier Madre
<carlosbo...@gmail.com> wrote:
> Hi,
> I've got something like f(Beta,eps)=g(
> integral(h(Beta,eps,tita1,tita2,...titaN)dtita1*dtita2*....*dtitaN), Beta,
> eps) where g and h are functions, and the integral of h showed as a variable
> of g has no primitive.
> I want to plot f, with Beta and eps in de x and y axes. Now I am just
> trying to plot with eps fixed and Beta in the X axis, but I get
> "PolynomialDivisionFailed" when it comes to plotting, I have followed "the
> calls" trying to ferret out how to fix it, but I am kind of a beginner and
> I'm a bit overwhelmed. I would really appreciate your help.
>
> Thanks in advance! Have a nice day:)
>
>
> This is my code:
>
> from numpy import *
> from sympy import *

I don't know if this is the source of your issue, but you should avoid
doing this. If you must mix SymPy and NumPy functions, use

import sympy as sym
import numpy as np

and use the functions like sym.cos or np.cos.

This will import all NumPy functions and then import all SymPy
functions, which overwrites the ones that have the same names as the
NumPy functions that were just imported. This is bad because SymPy
functions and NumPy functions do not mix very well. SymPy functions do
not work on NumPy arrays, and NumPy functions do not work on SymPy
expressions.

In general, it's a good idea to avoid import *, unless you are working
interactively (and even then, never do import * from more than one
library that may share the same function names).

Aaron Meurer

[Carlos Bouthelier Madre]

Thanks Aaron, didn't know that. I will fix it