$L$=$1 mH$
$E$=$2 V$
$R$=$1.0 \Omega$
$i(t)=i(t-\varDelta T)+\frac{1}{L}\int\limits_{t-\varDelta T}^{t} (-Ri(t)+E)\varDelta \tau$
$I(k) = (1-T_s\frac{R}{L}) I(k-1) + \frac{T_s}{L}E$
$I(k) = \frac{1}{(1+T_s\frac{R}{L})} I(k-1) + \frac{\frac{T_s}{L}}{(1+T_s\frac{R}{L})}E$
$I(k) = \frac{1-T_s\frac{R}{2L}}{(1+T_s\frac{R}{2L})} I(k-1) + \frac{\frac{T_s}{L}}{(1+T_s\frac{R}{2L})}E$
import numpy as np
import matplotlib.pyplot as plt
import math as math
np.set_printoptions(sign=' ')
# Circuit parameters:
E = 2.0
L = 1.0e-3
R = 1.0
# Total simulation time
T_total = 0.01
# 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))
# Solution vectors for current I
I_analytical = np.zeros((1, npoint))
I_ef = np.zeros((1, npoint))
I_eb = np.zeros((1, npoint))
I_tr = np.zeros((1, npoint))
# Time loop
# Initial condition
for i in np.arange(1,npoint):
# Analytical solution
I_analytical[0,i] = E/R*(1-math.exp(-R/L*i*Ts))
# Euler Forward used for discretization
I_ef[0,i] = (1-Ts*R/L) * I_ef[0,i-1] + Ts/L*E
# Euler Backward used for discretization
I_eb[0,i] = 1/(1+Ts*R/L) * I_eb[0,i-1] + Ts/L/(1+Ts*R/L)*E
# Trapezoidal Rule used for discretization
I_tr[0,i] = (1-Ts*R/(2*L))/(1+Ts*R/(2*L)) * I_tr[0,i-1] + Ts/L/(1+Ts*R/(2*L))*E
# Plots
# Time vector
t = np.arange(0, npoint)*Ts
plt.figure(figsize=(8,6))
plt.plot(t,I_analytical[0,:], linewidth=2)
plt.plot(t,I_ef[0,:],'--', linewidth=2)
plt.plot(t,I_eb[0,:],'--', linewidth=2)
plt.plot(t,I_tr[0,:],'--', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.legend(['Analtyical solution', 'Euler Forward', 'Euler Backward', 'Trapezoidal Rule'], prop={'size': 14})
plt.show()