Bottom Page

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.
Quote

Top Page

Possibly Related Threads...
Thread Author Replies Views Last Post
  Converting complex codebase from 2.7 to 3.5 rlgoodman 3 140 Apr-18-2019, 09:17 PM
Last Post: francisco_neves2020
  AssertionError: View function mapping is overwriting an existing endpoint function Zhavi221 7 285 Apr-17-2019, 01:07 PM
Last Post: Zhavi221
  problem with complex number jodrickcolina 1 90 Apr-13-2019, 06:59 PM
Last Post: Yoriz
  parsing complex text file anna 1 130 Apr-10-2019, 09:54 PM
Last Post: Larz60+
  Complex Integer saumyagoel 2 233 Feb-20-2019, 12:01 AM
Last Post: scidam
  Complex Long codes ayaz786amd 5 343 Dec-21-2018, 02:37 PM
Last Post: ichabod801
  Slicing a complex dictionary and list Drone4four 2 388 Aug-16-2018, 02:27 PM
Last Post: ichabod801
  How to detect semantic errors with numbers more complex than integers? RedSkeleton007 1 365 Aug-01-2018, 09:27 AM
Last Post: Larz60+
  Decypher complex (to me) statement PreservedKillich 5 548 Jul-31-2018, 06:13 PM
Last Post: PreservedKillich
  Run function in parallel but inherite dynamic data from the previous function atizva 4 567 Jul-11-2018, 06:39 AM
Last Post: volcano63

Forum Jump:


Users browsing this thread: 1 Guest(s)