Computational Science Asked by Tomás Lopes on February 9, 2021
I am trying to perform the following integral
$$int_{0}^{2pi}int_{0}^{+infty} frac{r’left(e^{-r’^2/2sigma^2}right)left(r-r’cos(theta-theta’)right)}{r^2+r’^2-2rr’cos(theta-theta’)}dr’dθ’$$
Using Gauss-Hermite for $r$ and Simpson 1/3 rule for $theta$ with no success. I can’t find my mistake but the output should look like Fig. 2. This was my code (sorry for my bad formatting, this is my first time uploading here).
$sigma$ should be assumed as 1.
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as ss
def rt(d, r, theta ,sig):
return r*(d-r*np.cos(theta))*np.exp(-r**2/(2*sig**2))/(d**2+r**2-2*d*r*np.cos(theta))
def intheta1(d, r, b, sig, N):
h = b/N
I = rt(d,r,0,sig) + rt(d,r,b,sig)
for i in range(1, N, 2):
I += 4*rt(d, r, i*h, sig)
for j in range(2, N, 2):
I += 2*rt(d, r, j*h, sig)
return I*h/3
def intr1(d, b, sig, N, M):
x, w = ss.roots_hermitenorm(N)
s = 0
for k in range(N):
s += intheta1(d, x[k], b, sig, M)*w[k]
return s/2
ps = np.linspace(0, 5, 1000)
qs = intr1(xs, 2*np.pi, 1, 1000, 90)
plt.plot(ps, qs)
I have been studying this type of numerical integration and I believe I understood my mistake. First of all I am using Gauss-Hermite which work with limits ${-infty}$ to ${infty}$ so using the fact that this function is even makes it so that to integrate from $0$ to ${infty}$ I have to use np.abs()
of my integration variable. Also, using Gauss-Hermite makes it so that I have to remove the exponential function. In this case I am using roots_hermitenorm()
so I had to find a way to remove $exp(-r^2/2)$ from the expression.
I got to these answer which is currently working flawlessly.
python
def integral_theta(r, rline, theta, sigma):
rline = np.abs(rline)
return np.exp((-(rline)**2*(1-sigma**2))/2/sigma**2) * rline * (r - rline*np.cos(theta))/(r**2 + rline**2 - 2*r*rline*np.cos(theta))
def i_theta(r, rline, sigma):
a, b = 0, 2*np.pi
N = 100
h = (b-a)/N
s_odd = 0
for k in range(1,N,2):
s_odd += integral_theta(r, rline, a+k*h, sigma)
s_even = 0
for j in range(2, N-1,2):
s_even += integral_theta(r, rline, a+j*h, sigma)
return h/3*(integral_theta(r, rline, a, sigma) + integral_theta(r, rline, b, sigma) + 4*s_odd + 2*s_even)
def i_r(r, sigma):
M = 1000
x, w = ss.roots_hermitenorm(M)
s = 0
for h in range(M):
s += i_theta(r, x[h], sigma)*w[h]
return s/2/np.pi/sigma**2
Answered by Tomás Lopes on February 9, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP