$R_1$=$1 \Omega$, $R_2$=$2 \Omega$
$L_1$=$20 mH$, $L_2$=$100 mH$, $L_3$=$50 mH$
$V_{in}(t)$=$10 V sin(\omega t)$
import numpy as np
import ipywidgets as widget
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
np.set_printoptions(sign=' ')
# Circuit parameters
R1 = 1
R2 = 2
L1 = 20e-3
L2 = 100e-3
L3 = 50e-3
Emax = 10
G1 = 1/R1
G2 = 1/R2
#default start values for Ts and T_total
Ts=1e-7
T_total=0.05
#function to perform the euler backward integration
def int_EB(Ts, T_total, npoint):
GL1_back = Ts/(L1)
GL2_back = Ts/(L2)
GL3_back = Ts/(L3)
# Conductance matrix
Gn = np.array([ [G1+GL1_back, -GL1_back, 0],
[-GL1_back, GL1_back+GL2_back+GL3_back, -GL3_back],
[0, -GL3_back, GL3_back+G2]])
# Node voltage vector
vn_back = np.zeros((3,npoint))
# Current vector in the 3 inductances
in_back = np.zeros((3,npoint))
# Intial conditions
# Voltage source at E(t=0)=0;
# i_L(t=0)=0 -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
in_back[:,0] = np.zeros(3)
vn_back[:,0] = np.zeros(3)
# Enter Loop
# tic %Start a timer
for i in np.arange(1,npoint):
# Update source vector
AL1_back = in_back[0,i-1]
AL2_back = in_back[1,i-1]
AL3_back = in_back[2,i-1]
E = Emax*np.sin(2*np.pi*60*i*Ts)
Jn = np.array([E/R1-AL1_back, AL1_back-AL2_back-AL3_back, AL3_back])
# Matrix inversion and solution of the equation G*e=J;
vn_back[:,i] = np.linalg.solve(Gn, Jn)
# post step
in_back[0,i] = AL1_back+GL1_back*(vn_back[0,i]-vn_back[1,i])
in_back[1,i] = AL2_back+GL2_back*(vn_back[1,i])
in_back[2,i] = AL3_back+GL3_back*(vn_back[1,i]-vn_back[2,i])
# toc %stop the timer
return in_back
#function to perform the euler forward integration
def int_EF(Ts, T_total, npoint):
# Large resistor added in parallel with inductors (3 current sources coincide at a node)
Gadd = 1e-4
# Conductance matrix
Gn = np.array( [[G1+Gadd, -Gadd, 0],
[-Gadd, 3*Gadd, -Gadd],
[ 0, -Gadd, Gadd+G2]])
# Node voltage vector
vn_forw = np.zeros((3,npoint))
# Current vector in the 3 inductances
in_forw = np.zeros((3,npoint))
# Intial conditions
# Voltage source at E(t=0)=0;
# i_L(t=0)=0 -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
in_forw[:,0] = np.zeros(3)
vn_forw[:,0] = np.zeros(3)
#Enter Loop
# tic % Start a timer
for i in np.arange(1,npoint):
# Update source vector
AL1_forw = in_forw[0,i-1]+(vn_forw[0,i-1]-vn_forw[1,i-1])*Ts/L1
AL2_forw = in_forw[1,i-1]+(vn_forw[1,i-1])*Ts/L2
AL3_forw = in_forw[2,i-1]+(vn_forw[1,i-1]-vn_forw[2,i-1])*Ts/L3
E = Emax*np.sin(2*np.pi*60*i*Ts)
Jn = np.array([E/R1-AL1_forw, AL1_forw-AL2_forw-AL3_forw, AL3_forw])
# Matrix inversion and solution of the equation G*e=J
vn_forw[:,i] = np.linalg.solve(Gn, Jn)
# post step
in_forw[0,i]= AL1_forw
in_forw[1,i]= AL2_forw
in_forw[2,i]= AL3_forw
return in_forw
#function to perform the trapezoidal integration
def int_TR(Ts, T_total, npoint):
GL1_trap = Ts/(2*L1)
GL2_trap = Ts/(2*L2)
GL3_trap = Ts/(2*L3)
# Conductance matrix
Gn = np.array( [[G1+GL1_trap, -GL1_trap, 0],
[-GL1_trap, GL1_trap+GL2_trap+GL3_trap, -GL3_trap],
[0, -GL3_trap, GL3_trap+G2]])
# Node voltage vector
vn_trap = np.zeros((3,npoint))
# Current vector in the 3 inductances
in_trap = np.zeros((3,npoint))
# Intial conditions
# Voltage source at E(t=0)=0
# i_L(t=0)=0 -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
in_trap[:,0] = np.zeros(3)
vn_trap[:,0] = np.zeros(3)
#Enter Loop
#tic
for i in np.arange(1,npoint):
# Update source vector
AL1 = in_trap[0,i-1]+(vn_trap[0,i-1]-vn_trap[1,i-1])*Ts/(2*L1)
AL2 = in_trap[1,i-1]+(vn_trap[1,i-1])*Ts/(2*L2)
AL3 = in_trap[2,i-1]+(vn_trap[1,i-1]-vn_trap[2,i-1])*Ts/(2*L3)
E = Emax*np.sin(2*np.pi*60*i*Ts)
Jn = np.array([E/R1-AL1, AL1-AL2-AL3, AL3])
# Matrix inversion and solution of the equation G*e=J
vn_trap[:,i] = np.linalg.solve(Gn, Jn)
# post step
in_trap[0,i]= AL1+GL1_trap*(vn_trap[0,i]-vn_trap[1,i])
in_trap[1,i]= AL2+GL2_trap*(vn_trap[1,i])
in_trap[2,i]= AL3+GL3_trap*(vn_trap[1,i]-vn_trap[2,i])
#toc
return in_trap
#Plot function
def plot_fun(time_step, npoint):
#Number of subplots (1*3) and size of the figure where the subplots will appear
fig, ax = plt.subplots(1,3, figsize=(16, 5))
#Calculate time values
t = np.arange(0, npoint)*Ts
#Transform time values to milliseconds for better representation
t= np.multiply(t, 10e+3)
#Plot results E. Backward
ax[0].set_title('E. Backward')
ax[0].set_xlim([0, (npoint-1)*Ts*10e+3])
ax[0].set_xlabel('Time [ms]')
ax[0].set_ylabel('Current [A]')
ax[0].plot(t,in_back[0,:], 'r', t, in_back[2,:], 'r:')
leg = ax[0].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
#Plot results E. Forward
ax[1].set_title('E. Forward')
#Set axis limite and transform to milliseconds
ax[1].set_xlim([0, (npoint-1)*Ts*10e+3])
ax[1].set_xlabel('Time [ms]')
ax[1].set_ylabel('Current [A]')
ax[1].plot(t, in_forw[0,:] ,'g', t, in_forw[2,:], 'g:', linewidth=2)
leg = ax[1].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
#Plot results Trapezoidal
ax[2].set_title('Trapezoidal')
ax[2].set_xlim([0, (npoint-1)*Ts*10e+3])
ax[2].set_xlabel('Time [ms]')
ax[2].set_ylabel('Current [A]')
ax[2].plot(t,in_trap[0,:], 'b', t, in_trap[2,:], 'b:', linewidth=2)
leg = ax[2].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
#Slider function to forward the selected values to the simulation functions
def slider_fun(time_step):
global in_back , in_forw , in_trap, npoint
if time_step==0 or T_total==0:
print('\n \n \n' + 'Simulation time step is set to: ' + str(time_step) + ' seconds')
print('Total simulation time: ' + str(T_total) + ' seconds' + '\n \n')
else:
print('\n \n \n' + 'Simulation time step is set to: ' + str(time_step) + ' seconds')
print('Total simulation time: ' + str(T_total) + ' seconds' + '\n \n')
npoint = int(np.round(T_total/time_step))
in_back=int_EB(time_step, T_total, npoint)
in_forw=int_EF(time_step, T_total, npoint)
in_trap=int_TR(time_step, T_total, npoint)
plot_fun(time_step, npoint)
return time_step, T_total
#Values range for the simulation steps
values=[round(i*10**-7, 7) for i in range(101)]
#It is important to set the argument continous_update to 'false' so that the slider_fun is called after the user finishes dragging the slider.
#Otherwise the function will be continuously called for several values in the dragging interval and the expected results will take alot of time to appear
slider= widget.interactive(slider_fun, time_step=widget.SelectionSlider(description="Time step", continuous_update=False, options=[("%g"%i,i) for i in values], layout=widget.Layout( width='50%', height='50px',)));
display(slider)
$7,6 \cdot 10^{-6} s < T_{critical} $