$R$=$1 \Omega$
$R_1 = R_2 = R_3 = R$
$L$=$0.5 mH$
$C$=$0.5 mF$
$C_1 = C_2 = C_3 = C$
$I_3$=$1 A$
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(sign=' ')
# Circuit parameters:
R = 1.0
L = 0.5e-3
C = 0.5e-3
I3 = 1.0
# Incidence matrix
M = np.array([ [1, 1, 0],
[-1, 0, 1],
[0, -1, -1]])
# Number of branches Nb and number of nodes Nn:
(Nn, Nb) = np.shape(M)
# Branches
R_b = np.array([ [R, 0.5*R, R]])
L_b = np.array([ [0.5*L, L, 0.5*L]])
# Nodes
R_n = np.array([ [R, R, R]])
C_n = np.array([ [C, C, C]])
# Total simulation time
T_total = 0.008
# Simulaiton time step
Ts = 0.1e-3
# Number of simulation time steps
npoint = int(np.round(T_total/Ts))
print('Total simulation time: ' + str(T_total))
print('\nSimulation time step: ' + str(Ts))
# Stability condition
print('\nStability condition:\n Ts < sqrt(Lk*Ck)')
# Calculate maximum allowed time step based on circuit parameters
Ts_max = np.sqrt(np.amin([ [L_b[0,0]*C_n[0,0]], [L_b[0,0]*C_n[0,1]], [L_b[0,1]*C_n[0,0]], [L_b[0,1]*C_n[0,2]], [L_b[0,2]*C_n[0,1]], [L_b[0,2]*C_n[0,2]]]))
print(str(Ts) + ' < ' + str(Ts_max))
# Current source vectors for all nodes
H = np.zeros((npoint, Nn))
H[:,2] = I3*np.ones((npoint))
# Voltage source vectors for all branches
E = np.zeros((npoint, Nb))
# Create matrix to store solution of branch currents in every simulation steps
I = np.zeros((npoint, Nb))
# Create matrix to store solution of node voltages in every simulation steps
V = np.zeros((npoint, Nn))
# Initial conditions
# Branch inductances: current = 0
I[0,:] = np.zeros((1, Nb))
# Node capacitance: voltage = 0
# Calculate V(Ts/2) based on a given initial conditions V(0)
V_0 = np.zeros((1, Nn))
for j in np.arange(0, Nn):
V[0,j] = (C_n[0,j]*V_0[0,j]/Ts - np.matmul(M[j,:], np.transpose(I[0,:])) + H[0,j]) / (C_n[0,j]/Ts + 1/R_n[0,j])
# Time loop
for i in np.arange(1,npoint):
# Branch loop - calculation of I(i*Ts) based on I((i-1)*Ts) and V((i-1/2)*Ts)
for j in np.arange(0, Nb):
I[i,j] = I[i-1,j] + Ts/L_b[0,j] * (np.matmul(np.transpose(M[:,j]), np.transpose(V[i-1,:])) - R_b[0,j] * I[i-1,j] + E[i-1,j])
# Node loop - calculation of V((i+1/2)*Ts) based on I(i*Ts) and V((i-1/2)*Ts)
for j in np.arange(0, Nn):
V[i,j] = (C_n[0,j]*V[i-1,j]/Ts - np.matmul(M[j,:], np.transpose(I[i,:])) + H[i,j]) / (C_n[0,j]/Ts + 1/R_n[0,j])
# Time vectors
# for branch currents
t_I = np.arange(0, npoint)*Ts
# for node voltages
t_V = (np.arange(0, npoint) + 0.5)*Ts
# Plots
plt.figure(figsize=(8,6))
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.title('Branch currents')
plt.plot(t_I,I[:,0], t_I, I[:,1], t_I,I[:,2], linewidth=2)
plt.xlim([0, (npoint-1)*Ts])
plt.legend(['Current $i_{a}$(t)', 'Current $i_{b}$(t)', 'Current $i_{c}$(t)'], prop={'size': 14})
plt.show()
# Plots
plt.figure(figsize=(8,6))
plt.xlabel('Time [s]')
plt.ylabel('Voltage [V]')
plt.title('Node voltages')
plt.plot(t_V,V[:,0], t_V, V[:,1], t_V,V[:,2], linewidth=2)
plt.xlim([0, (npoint-1)*Ts])
plt.legend(['Voltage $v_{1}$(t)', 'Voltage $v_{2}$(t)', 'Voltage $v_{3}$(t)'], prop={'size': 14})
plt.show()