-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathBessel.py
96 lines (82 loc) · 2.25 KB
/
Bessel.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#Feito por Hevenicio Silva
# Importação de biblitecas
from IPython.display import display, Markdown, Latex
import matplotlib.pyplot as plt
import pandas as pd
import sympy as sym
import numpy as np
# Entrada de dados
dados = pd.read_csv("dados.csv")
df = pd.DataFrame(dados)
# Função Fatorial
def fatorial(n):
if(n == 0):
return 1
else:
return n*fatorial(n - 1)
# Formatação dos dados impressos - Latex
def printmd(string):
display(Markdown(string))
#display(Latex(string))
f = sym.symbols('$_{0}$')
g = sym.symbols('$_{1}$')
h = sym.symbols('$_{2}$')
# Definindo as Funções de Bessel
def J0(x):
k = 0
soma = 0
termo = 1
while(1.0e-4 <= termo):
termo = pow(x/2, (2*k))/(fatorial(k)*fatorial(k))
soma += pow(-1, k)*termo
k += 1
return soma
def J1(x):
k = 0
soma = 0
termo = 1
while(1.0e-4 <= termo):
termo = pow(x/2, (2*k + 1))/(fatorial(k)*fatorial(k + 1))
soma += pow(-1, k)*termo
k += 1
return soma
def J2(x):
k = 0
soma = 0
termo = 1
while(1.0e-4 <= termo):
termo = pow(x/2, (2*k + 2))/(fatorial(k)*fatorial(k + 2))
soma += pow(-1, k)*termo
k += 1
return soma
# Amplitude da função de Bessel
y = np.sqrt(2/(np.pi*df['entrada']))
z = -np.sqrt(2/(np.pi*df['entrada']))
# Resultados obtidos
l =[]
m = []
n = []
o = []
p = [l, m, n, o]
for i in df['entrada']:
l.append(i)
m.append(J0(i))
n.append(J1(i))
o.append(J2(i))
#printmd('**$J${}({}): {: 0.5f} $|$ $J${}({}): {: 0.5f} $|$ $J${}({}): {: 0.5f}**'.format(f, i, J0(i), g, i, J1(i), h, i, J2(i)))
a = pd.DataFrame(p)
a = a.rename(index={0:"x", 1:"J0(x)", 2: "J1(x)", 3:"J2(x)"}).T
print(a)
# Plotando os resultados obtidos
plt.title('Funções de Bessel - J$_{0}(x)$, J$_{1}(x)$ e J$_{2}(x)$')
plt.plot(df['entrada'], m, 'g', label = 'J$_{0}(x)$')
plt.plot(df['entrada'], n, 'limegreen', label = 'J$_{1}(x)$')
plt.plot(df['entrada'], o, 'lawngreen', label = 'J$_{2}(x)$')
plt.plot(df['entrada'], y, 'c--', label = '$\sqrt{2 / \pi x}$')
plt.plot(df['entrada'], z, 'b--', label = '$ - \sqrt{2 / \pi x}$')
plt.legend(framealpha=1, frameon=True)
plt.grid(color = 'black')
plt.style.use('seaborn')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()