#import packages needed
import numpy as np
import scipy as sp
import scipy.integrate
import matplotlib.pyplot as plt
# Solve ODE numerically
# start with Rv = 40 kiloohms to observe the region where we expect a single periodic solution
def dx_dt(t, x):
R = 47
Rv = 40
V0 = 0.250
R0 = 157
R2 = 6
R1 = 1
x1 = x[0]
x2 = x[1]
x3 = x[2]
return (x2, x3, (-(R / Rv) * x3) - x2 - ((R2 / R1) * min(x1, 0)) - ((R / R0) * V0))
ts = np.linspace(100, 200, 1000)
res = sp.integrate.solve_ivp(dx_dt, [0, 200], y0 = [0.25, 0, 0], t_eval = ts, rtol = 1e-9)
print("Solution to Third Order Differential Equation for R_v = 40 kiloohms")
plt.figure(figsize = (9, 6))
plt.plot(ts, res.y[0], '-', color = 'black', label = 'x1')
# plt.plot(ts, res.y[1], '-', color = 'red', label = 'x2')
# plt.plot(ts, res.y[2], '-', color = 'blue', label = 'x3')
plt.title('Fig.1 \nTime Series \nX')
plt.ylabel('Voltage (volts)')
plt.xlabel('Time (s)')
# plt.legend()
plt.show()
# determine the maximum amplitude
max_value = np.amax(res.y[0])
min_value = np.amin(res.y[0])
max_amplitude = (abs(max_value) + abs(min_value)) / 2
print(f'Maximum amplitude for the solution at Rv = 40 kiloohms: {max_amplitude:.2f}V')
# try Rv = 60 kiloohms to observe the region where we expect it to bifurcate and have two periodic solutions
def dx_dt(t, x):
R = 47
Rv = 60
V0 = 0.250
R0 = 157
R2 = 6
R1 = 1
x1 = x[0]
x2 = x[1]
x3 = x[2]
return (x2, x3, (-(R / Rv) * x3) - x2 - ((R2 / R1) * min(x1, 0)) - ((R / R0) * V0))
ts = np.linspace(100, 200, 1000)
res = sp.integrate.solve_ivp(dx_dt, [0, 200], y0 = [0.25, 0, 0], t_eval = ts, rtol = 1e-9)
print("Solution to Third Order Differential Equation for R_v = 60 kiloohms")
plt.figure(figsize = (9, 6))
plt.plot(ts, res.y[0], '-', color = 'black', label = 'x1')
# plt.plot(ts, res.y[1], '-', color = 'red', label = 'x2')
# plt.plot(ts, res.y[2], '-', color = 'blue', label = 'x3')
plt.title('Fig.2 \nTime Series \nX')
plt.ylabel('Voltage (volts)')
plt.xlabel('Time (s)')
# plt.legend()
plt.show()
# determine the maximum amplitude
max_value = np.amax(res.y[0])
min_value = np.amin(res.y[0])
max_amplitude = (abs(max_value) + abs(min_value)) / 2
print(f'Maximum amplitude for the solution at Rv = 60 kiloohms: {max_amplitude:.2f}V')
# try Rv = 75 kiloohms to observe the region where we expect chaotic behavior (inifinite solutions)
def dx_dt(t, x):
R = 47
Rv = 75
V0 = 0.250
R0 = 157
R2 = 6
R1 = 1
x1 = x[0]
x2 = x[1]
x3 = x[2]
return (x2, x3, (-(R / Rv) * x3) - x2 - ((R2 / R1) * min(x1, 0)) - ((R / R0) * V0))
ts = np.linspace(100, 200, 1000)
res = sp.integrate.solve_ivp(dx_dt, [0, 200], y0 = [0.25, 0, 0], t_eval = ts, rtol = 1e-9)
print("Solution to Third Order Differential Equation for R_v = 75 kiloohms")
plt.figure(figsize = (9, 6))
plt.plot(ts, res.y[0], '-', color = 'black', label = 'x1')
# plt.plot(ts, res.y[1], '-', color = 'red', label = 'x2')
# plt.plot(ts, res.y[2], '-', color = 'blue', label = 'x3')
plt.title('Fig.3 \nTime Series \nX')
plt.ylabel('Voltage (volts)')
plt.xlabel('Time (s)')
# plt.legend()
plt.show()
# determine the maximum amplitude
max_value = np.amax(res.y[0])
min_value = np.amin(res.y[0])
max_amplitude = (abs(max_value) + abs(min_value)) / 2
print(f'Maximum amplitude for the solution at Rv = 75 kiloohms: {max_amplitude:.2f}V')