-
Notifications
You must be signed in to change notification settings - Fork 0
/
EDO - Oscilador Harmônico.py
131 lines (113 loc) · 4.51 KB
/
EDO - Oscilador Harmônico.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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
# made by Leonardo Saads
import matplotlib.pyplot as plt
import math
import numpy as np
# Classe dos Erros
class ValorNegativoError(Exception):
def __init__(self):
super().__init__("Erro de Valor Menor ou Igual a Zero")
class Oscilador:
"""
Implementa a classe do Oscilador Harmônico Amortecido
"""
def __init__(self, m: float, u: float, k: float, x_inicial: float, v_inicial):
"""
Adquire os parâmetros da classe.
:param m: massa
:param u: constante positiva
:param k: constante elástica
:param x_inicial: espaço inicial
:param v_inicial: valocidade inicial
"""
if m <= 0 or u <= 0 or k <= 0:
raise ValorNegativoError
else:
self.m = m
self.u = u
self.k = k
self.space = x_inicial
self.velocity = v_inicial
def __repr__(self):
"""
Representação.
:return: string
"""
return f'------- SISTEMA OSCILADOR -------\n' \
f'Valor da Massa: {self.m}\n' \
f'Valor do U: {self.u}\n' \
f'Valor do K: {self.k}\n' \
f'Condições Iniciais do Sistema' \
f'f(0) = {self.space}\n' \
f'f`(0) = {self.velocity}\n'
def lambdas(self):
"""
Implelemta o cálculo das raízes dos polinômio característico.
:return: lambda1, lambda2
"""
v = (self.u/self.m)/2
w = (self.k/self.m)**(1/2)
lambda1 = (v**2 - w**2)**(1/2) - v
lambda2 = - (v**2 - w**2)**(1/2) - v
return lambda1, lambda2
def funcao(self):
"""
determina a funcao a ser utilizada
:return: funcao
"""
raiz1, raiz2 = Oscilador.lambdas(self)
if raiz1 == raiz2:
return lambda t: self.space*math.e**(raiz1*t) + (self.velocity - self.space*raiz1)*t*math.e**(raiz2*t)
if raiz1 != raiz2 and type(raiz2) == float:
coef_2 = (self.velocity - self.space * raiz1) / (raiz2 - raiz1)
coef_1 = (self.velocity - self.space * raiz2) / (raiz1 - raiz2)
return lambda t: coef_1*math.e**(raiz1*t) + coef_2*math.e**(raiz2*t)
else:
return lambda t: self.space*math.e**(raiz1.real*t)*math.cos(raiz1.imag*t) + \
((self.velocity+raiz2.real*self.space)/raiz2.imag)*math.e**(raiz2.real*t) *\
math.sin(raiz2.imag*t)
def aproximacao(self, h: float, final: float):
"""
A aproximação de uma EDO com base no método de euler.
:return: vetor
"""
# z'[n+1] = z[n] + h*w[i-1]
# w'[n+1] = w[n] + h*(-2v*w[n] - w**2 z[n])
# y(0) = z(0) & y'(0) = w(0)
# defininindo z e w
v = (self.u / self.m) / 2
g = (self.k / self.m) ** (1 / 2)
z = np.zeros(len(np.arange(0, final, h))) # um vetor de zeros
w = np.zeros(len(np.arange(0, final, h))) # um vetor de zeros
z[0] = self.space # valor inicial
w[0] = self.velocity # valor incial
for i in range(1, int(final/h)):
z[i] = z[i-1] + h*w[i-1]
w[i] = w[i-1] + h*((-2*v)*z[i-1] - (g**2)*w[i-1])
return z, w
def grafico(self, largura: int):
"""
implementa o gráfico
:largura: a largura do gráfico.
:return: gráfico
"""
# parte de vetores - aproximações e euler
aproximacoes = Oscilador.aproximacao(self, 0.05, largura)[0] # DETERMINE AQUI A PRECISÃO DA APROXIMAÇÃO
x_aproc = np.linspace(0, largura, len(aproximacoes))
# Partes gráficas
plt.plot(x_aproc, aproximacoes, label='Euler method aproximation')
plt.xlabel('Tempo t')
plt.ylabel('Space x(t)')
plt.legend()
plt.grid(True)
plt.show()
# exemplo de chamada da classe
if __name__ == '__main__':
mass = float(input("Detemine o valor da Massa: \n"))
uconst = float(input("Detemine o valor da U: \n"))
kconst = float(input("Detemine o valor do K: \n"))
y0 = float(input("determine as condições iniciais do Sistema - y(0): \n"))
y10 = float(input("determine as condições iniciais do Sistema - y'(0): \n"))
oscilador_teste = Oscilador(mass, uconst, kconst, y0, y10)
print('Raízes do polinômio caracterítico')
print(oscilador_teste.lambdas())
oscilador_teste.grafico(30)