SymPy: использование решающего набора для решения интеграла - PullRequest
Новая
       15

SymPy: использование решающего набора для решения интеграла

0 голосов
/ 14 марта 2019

Я пытаюсь использовать solveset для решения отношения интегралов, при этом решением является разрешенная область комплексного числа.

Я написал: 

import numpy as np
import sympy
import cmath
import scipy.integrate as integrate
from sympy import symbols, Eq, I, solveset

m1 = 5.28
mlep = 0.11
mtau = 1.78 
m2 = 1.86
q2_max = (m1-m2)**2


def p2(q2):
    return np.sqrt(((m1-m2)**2 - q2)*((m1+m2)**2 - q2))

def s(q2):
    return q2/m1**2

def fplus(q2):
    return 0.79/(1-0.75*s(q2) + 0.039*s(q2)**2)
def fminus(q2):
    return -0.36/(1-0.77*s(q2)-0.098*(s(q2)**2))
def fs(q2):
    return 0.80/(1-0.22*s(q2) - 0.098*s(q2)**2)
def ft(q2):
    return 0.77/(1-0.76*s(q2) + 0.043*s(q2)**2)

def h0(q2):
    return (2*m1*p2(q2)*fplus(q2))/(np.sqrt(q2))
def ht(q2):
    return (1/np.sqrt(q2))*((m1**2 - m2**2)*fplus(q2)+q2*fminus(q2))
def hsp(q2):
    return (m1+m2)*fs(q2)
def h_t(q2):
    return (2*m1*p2(q2)*ft(q2))/(m1+m2)

def lep(q2):
    return p2(q2)*q2*((1-((mlep**2)/q2))**2)*(np.absolute(h0(q2))**2)

lep_int = integrate.quad(lep, mlep**2, q2_max)[0]

def delta_tau(q2):
    return mtau**2/(2*q2)

def tau_np(q2,v_l,v_r,s_l,s_r,t_l):
    return p2(q2)*q2*((1-((mtau**2)/q2))**2)*((abs(1+v_l+v_r)**2)*(((1+delta_tau(q2))*(np.absolute(h0(q2))**2))+(3*delta_tau(q2)*(np.absolute(ht(q2))**2)))+((3/2)*(np.absolute(s_l+s_r)**2)*(np.absolute(hsp(q2))**2))+(8*np.absolute(t_l)**2*(1+4*delta_tau(q2))*np.absolute(h_t(q2))**2)+(3*np.sqrt(2*delta_tau(q2))*(s_l+s_r).real*hsp(q2)*ht(q2))+(12*np.sqrt(2*delta_tau(q2))*t_l.real*h0(q2)*h_t(q2)))

def tau_int_np(v_l,v_r,s_l,s_r,t_l):
    return integrate.quad(tau_np, mtau**2, q2_max, args=(v_l,v_r,s_l,s_r,t_l))[0]

def r_d_np(v_l,v_r,s_l,s_r,t_l):
    return tau_int_np(v_l,v_r,s_l,s_r,t_l)/lep_int

x,y=symbols('x y', real=True)
v_l = x+I*y

expr = solveset(0.0<=r_d_np(v_l,0.0,0.0,0.0,0.0)<=1.0,(x,y))

Если я оцениваю r_d_np(0,0,0,0,0), это дает число около 0,115 - однако при попытке запустить это выдает ошибку quadpack.error: Supplied function does not return a valid float.

...