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
|
#!/usr/bin/env python
import matplotlib.pyplot as plt
import tkinter
from math import sin,cos,e
from sys import argv
def f(x):
return G*sin(A*x) + B*cos(B*x) + 2*e**(rer1 * x) * (reE * cos(imr1 * x) - imE * sin(imr1 * x))
for w in range(100,10000,100):
c1 = 2.5*10**(-6)
c2 = 1*10**(-6)
r1 = 200
rload = 1000
a = 1
b = c1*r1*rload + r1*c2*rload + (rload**2)*c2
c = rload
det = (b**2-4*a*c)**(1/2)
r1 = (-b-det)/(2*a)
r2 = (-b+det)/(2*a)
A = c1*c2*rload**2*w*(rload-w**2)/(b**2*w**2 + rload**2 - 2*rload*w**2 + w**4)
B = b*c1*c2*rload**2*w**3/(b**2*w**2 + rload**2 - 2*rload*w**2 + w**4)
C = - c1*c2*rload**2*w*(rload-w**2)/(b**2*w**2 + rload**2 - 2*rload*w**2 + w**4)
D = - b*c1*c2*rload**3*w/(b**2*w**2 + rload**2 - 2*rload*w**2 + w**4) # from https://www.wolframalpha.com/input/?i=solve+for+x1%2Cx2%2Cx3%2Cx4+in+%7B%7B1%2C0%2C1%2C0%7D%2C+%7Bb%2C1%2C0%2C1%7D%2C+%7BR%2C+b%2C+w%5E2%2C+0%7D%2C+%7B0%2C+R%2C+0%2C+w%5E2%7D%7D*%7Bx1%2Cx2%2Cx3%2Cx4%7D+%3D+%7B0%2C0%2Cw*c_1*c_2*R%5E2%2C0%7D and an insane partial fraction decomposition
E = (D-C*r1)/(r2-r1)
F = C - E # Another PFD of Cs-D/(as^2+bs+c)
G = A/w
# Final solution should be
# Gsin(A * theta) + Bcos(B * theta) + Ee^(r1 t) + Fe^(r2 t)
# But E, F, r1, and r2 are complex. Luckily, conjugates make it that
# = 2e^(Re(r1) t) ( Re(E)cos(Im(r1)t) - Im(E)sin(Im(r1)t) )
rer1 = r2.real
reE = F.real
imr1 = r2.imag
imE = F.imag # switched because positive ones were needed
print("%.2E * sin(%.2E t) + %.2E * cos(%.2E t) + 2e^(%.2f t) (%.2E cos(%.2f t) - %.2E sin(%.2f t))" % (G, A, B, B, rer1, reE, imr1, imE, imr1))
x = []
y = []
for i in range(1000):
x.append(i/100)
y.append(f(i/100))
plt.plot(x,y)
if len(argv) > 1 and argv[1] == 'img':
plt.savefig('plot.png')
else:
plt.show()
|
