Python Forum
Integration of a complex function having singularities using quad
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Integration of a complex function having singularities using quad
#1
Hi
In fact I want to visualize a function using python, but the function unfortunately includes singularities (one in it's end points and the other inside the boundary).
first we need to integrate the function f below from -1 to 1.
def f(u,t,x=l):
	q=u**2-x**2
	if q<0:
		q=-q
		v=complex(0,q)
	else:
		v=np.sqrt(u**2-x**2)	
	return(((1-u**2)/(u+v)**3)
then divide the integral of f by the integral of function f1.
def f1(u,t,x=l):
	q=u**2-x**2
	if q<0:
		q=-q
		v=complex(0,q)
	else:
		v=np.sqrt(u**2-x**2)		 
	return((u)*u*(1-((u-np.sqrt(v))/(u+np.sqrt(u**2+x**2)))**2))
Note:v is complex when q<0 to avoid negative number inside the square root.

We want to visualize the function that is the (integral of -1 to 1 of f divided by the integral of -1 t0 1 of f1 )times (4x**2) interms of x.(0<x<1)

first we need to use quad from SciPy in order to integrate both of the functions f1 and f. But as they are complex functions and quad cant integrate complex functions we need to define a function that does this by its own.

We can also notice that (\int_-1^1(f)/\int_-1^1(f1))times 4x^4 has singularities at -1 and 0. so we need to include options to let quad tackle these singularities.
so the code will be like:
import numpy as np
from scipy.integrate import quad
import scipy
import matplotlib.pyplot as plt
l=0  # variable representing x
ls=[] # list that represents the value of the function we want to visualize corresponding to certain value of l. 
lx=[] # list that includes all the values of l from 0 to 1.

for i in range(0,10000):
	ls.append([])
for i in range(0,10000):
	lx.append([])

def f(u,t,x=l):  #the third argument is a default argument
	q=u**2-x**2
	if q<0:
		q=-q
		v=complex(0,q)
	else:
		v=np.sqrt(u**2-x**2)	
	return((u-t)*(1-u**2)/(u+v)**3) #we multiply the function by (u-t) since quad by itself will multiply the function by 1/(u-wvar) to get the principle value of the integral 


def f1(u,t,x=l):
	q=u**2-x**2
	if q<0:
		q=-q
		v=complex(0,q)
	else:
		v=np.sqrt(u**2-x**2)		 
	return((u-t)*u*(1-((u-np.sqrt(v))/(u+np.sqrt(u**2+x**2)))**2))


def complex_quadrature(func,a,b):
 def real_func(c):
  return scipy.real(func(c,-1))
 def imag_func(c):
  return scipy.imag(func(c,-1))		
 real_integral=quad(real_func,a,b, weight='cauchy', wvar=-1)
 imag_integral=quad(imag_func,a,b,weight='cauchy',wvar=-1)  # we use this function to integrate from -1 to zero
 return(real_integral[0]+1j*imag_integral[0])
 
def complex_quadrature2(func,a,b):
 def real_func(c):
  return scipy.real(func(c,0))
 def imag_func(c):
  return scipy.imag(func(c,0))		
 real_integral=quad(real_func,a,b, weight='cauchy', wvar=0)    # we use this function to integrate from zero to 1
 imag_integral=quad(imag_func,a,b,weight='cauchy',wvar=0)
 return(real_integral[0]+1j*imag_integral[0])



s=(complex_quadrature(f,-1,0)+complex_quadrature2(f,0,1))/(complex_quadrature(f1,-1,0)+complex_quadrature2(f1,0,1))
  



for i in range(0,10000):
    ls[i]=4*(l**4)*s
    lx[i]=l                   #visualizing the function
    l=l+0.0001
	
	
plt.plot(ls,lx,color='b')
plt.show()	




the code will give the error:
Traceback (most recent call last):
  File "untitled11.py", line 48, in <module>
    s=(complex_quadrature(f,-1,0)+complex_quadrature2(f,0,1))/(complex_quadrature(f1,-1,0)+complex_quadrature2(f1,0,1)) 
  File "untitled11.py", line 27, in complex_quadrature
    real_integral=quad(real_func,a,b, weight='cauchy', wvar=-1)
  File "C:\Users\AGT-CR\AppData\Local\Programs\Python\Python37\lib\site-packages\scipy\integrate\quadpack.py", line 427, in quad
    raise ValueError(msg)
ValueError: Parameter 'wvar' must not equal integration limits 'a' or 'b'.
I tried also to make the endpoints be so close to 0 or 1 but their having their exact values.
s=(complex_quadrature(f,-0.9999,-0.0001)+complex_quadrature2(f,0.0001,1))/(complex_quadrature(f1,-0.9999.,-0.0001)+complex_quadrature2(f1,0.0001,1))  
but I will also get that error:
untitled11.py:21: RuntimeWarning: divide by zero encountered in double_scalars
  return((u-t)*(1-u**2)/(u+v)**3)
untitled11.py:27: IntegrationWarning: Extremely bad integrand behavior occurs at some points of the
  integration interval.
  real_integral=quad(real_func,a,b, weight='cauchy', wvar=-1)
untitled11.py:47: RuntimeWarning: divide by zero encountered in double_scalars
  return((u-t)*u*(1-((u-np.sqrt(v))/(u+np.sqrt(u**2+x**2)))**2))
C:\Users\AGT-CR\AppData\Local\Programs\Python\Python37\lib\site-packages\numpy\core\numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
So is there a way to handle such singularities? have I missed singular points other than 0 or -1?
Any help is really appreciated.
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  Active Directory integration dady 2 463 Oct-13-2023, 04:02 AM
Last Post: deanhystad
  Help with Integration Pandas excel - Python Gegemendes 5 1,723 Jun-05-2022, 09:46 PM
Last Post: Gegemendes
Photo Integration of apache spark and Kafka on eclipse pyspark aupres 1 3,701 Feb-27-2021, 08:38 AM
Last Post: Serafim
  Tableau Time Series Prediction using Python Integration tobimarsh43 0 1,887 Jul-24-2020, 10:38 AM
Last Post: tobimarsh43
  R-PYTHON INTEGRATION RELATED PROBLEM arnab93 0 1,418 Jun-05-2020, 02:07 PM
Last Post: arnab93
  STATA/Python Integration jprender 0 1,804 May-03-2020, 09:38 PM
Last Post: jprender
  phython language java integration jammytcs123123 1 2,265 Jul-04-2018, 03:13 AM
Last Post: Skaperen
  SQLite DB integration duplicate columns rachitmishra25 1 5,320 Oct-27-2017, 11:20 AM
Last Post: buran
  Trapezoidal integration method in python. Getting wrong results auting82 5 5,791 Oct-14-2017, 07:38 PM
Last Post: buran
  Sympy Integration Flexico 5 7,429 Dec-07-2016, 07:24 AM
Last Post: micseydel

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020