145 lines
6.5 KiB
Python
145 lines
6.5 KiB
Python
#!/usr/bin/python3
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from scipy.integrate import solve_ivp
|
|
# ==== Component values ====
|
|
R1 , C1 = 10e3 , 1e-6; tau1 = R1 * C1 # 0.010 s
|
|
R2 , C2 = 4.7e3 , 1e-6; tau2 = R2 * C2 # 0.0047 s
|
|
R4 , R5 = 1e3 , 4.7e3 ; G = 1 + R5 / R4 # 5.7
|
|
tau_min , tau_max = min ( tau1 , tau2 ) , max ( tau1 , tau2 )
|
|
# ==== System (cascade) ODE : x = [ v1 , v2 ] ====
|
|
def f_rhs (t , x , vin_func ) :
|
|
v1 , v2 = x
|
|
vin = vin_func (t)
|
|
dv1 = ( vin - v1 ) / tau1
|
|
dv2 = ( G * v1 - v2 ) / tau2
|
|
return [ dv1 , dv2 ]
|
|
|
|
# ==== Analytic impulse and step responses ====
|
|
def h (tt) :
|
|
tt = np.asarray ( tt , dtype = float )
|
|
out = np.zeros_like (tt)
|
|
pos = tt >= 0
|
|
# h (t) = G /( tau1 - tau2 ) * ( e ^{ - t / tau1 } - e ^{ - t / tau2 }) * u (t)
|
|
out [ pos ] = ( G /( tau1 - tau2 ) ) *( np.exp ( - tt [ pos ]/ tau1 ) - np.exp ( - tt [ pos ]/ tau2 ) )
|
|
return out
|
|
|
|
def s (tt) :
|
|
tt = np.asarray ( tt , dtype = float )
|
|
out = np.zeros_like (tt)
|
|
pos = tt >= 0
|
|
# s (t) = \ int_0 ^ t h (tau) d tu = G /( tau1 - tau2 ) [ tau1 (1 - e ^{ - t / tau1}) - tau2 (1 - e ^{ - t / tau2 }) ] u (t)
|
|
out [ pos ] = ( G /( tau1 - tau2 ) ) *( tau1 *(1.0 - np.exp ( - tt [ pos ]/ tau1 ) ) - tau2 *(1.0 - np.exp ( - tt [ pos ]/ tau2 ) ) )
|
|
return out
|
|
|
|
# ==== Two - stage time grid ( fine near pulse , coarse afterward ) ====
|
|
def make_time_grid ( T_pulse , t_end , fine_mult =3.0 , n_fine =400 , n_tail=2200) :
|
|
"""
|
|
Create a dense grid over [0 , fine_mult * T_pulse ] and a coarser
|
|
grid afterward to t_end .
|
|
"""
|
|
t_fine_end = min ( fine_mult * T_pulse , t_end )
|
|
t_fine = np.linspace (0.0 , t_fine_end , max (3 , n_fine ) )
|
|
t_tail = np.linspace ( t_fine_end , t_end , max (3 , n_tail ) )
|
|
# Ensure uniqueness and monotonicity ( remove duplicate junction point )
|
|
if t_tail.size and t_tail [0] == t_fine [ -1]:
|
|
t_tail = t_tail [1:]
|
|
return np.concatenate ([ t_fine , t_tail ]) if t_tail.size else t_fine
|
|
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
# Square pulse
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
A_rect , T_rect = 1.0 , 100e-6 # 1 V , 100 mu s
|
|
t_end_rect = 8.0* tau_max # long enough to see decay
|
|
|
|
# Input function ( square pulse on [0 , T_rect ) )
|
|
vrect = lambda t : A_rect if ( t >= 0.0 and t < T_rect ) else 0.0
|
|
|
|
# Two - stage t grid ( dense around the pulse )
|
|
t_rect = make_time_grid ( T_rect , t_end_rect , fine_mult =3.0 , n_fine=600 , n_tail =2400)
|
|
# Solve ODE ( v_in = square pulse )
|
|
sol_rect = solve_ivp ( f_rhs , ( t_rect [0] , t_rect [ -1]) , y0 =[0.0 , 0.0] ,
|
|
t_eval = t_rect , args =( vrect ,) ,
|
|
rtol =1e-8 , atol =1e-11 , method = 'RK45')
|
|
v2_rect = sol_rect.y [1]
|
|
# Area of the square pulse (exact)
|
|
S_rect = A_rect * T_rect
|
|
# Square pulse & overlay vs analytic impulse response
|
|
plt.figure ( figsize =(8 , 6) )
|
|
plt.subplot (2 ,1 ,1)
|
|
plt.plot ( t_rect , np.array ([ vrect (t) for t in t_rect ]) , 'k')
|
|
plt.grid (True) ; plt.xlabel ( 't (s) ') ; plt.ylabel ( 'V')
|
|
plt.title ( ' Section 1 - Square pulse input $v_{in}(t)$ ')
|
|
plt.subplot (2 ,1 ,2)
|
|
|
|
plt.plot ( t_rect , v2_rect / S_rect , label = r' $v_{\mathrm{ out }}/(AT)$ (ODE) ')
|
|
plt.plot ( t_rect , h(t_rect) , '--' , label = r' $h(t)$ (analytic)')
|
|
plt.grid (True) ; plt.xlabel ( 't (s) ') ; plt.ylabel ( 'V / ( V s ) ')
|
|
plt.title ( 'Area - normalized ODE response vs analytic impulse response')
|
|
plt.legend ( loc = 'best')
|
|
plt.tight_layout () ; plt.show ()
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
# Triangular pulse
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
A_tri , T_tri = 2.0 , 100e-6 # peak 2 V on [0 , T ] , area = 0.5* A * T = 1e-4 V s
|
|
t_end_tri = 8.0* tau_max
|
|
def vtri_scalar (t) :
|
|
""" Triangular pulse on [0 , T_tri ] with peak A_tri at T /2; 0 elsewhere."""
|
|
if t < 0.0 or t >= T_tri :
|
|
return 0.0
|
|
# symmetric triangle : v = (2* A / T ) * min (t , T - t )
|
|
|
|
return (2.0* A_tri / T_tri ) * min (t , T_tri - t )
|
|
|
|
# Dense grid around triangle ; coarser afterward
|
|
t_tri = make_time_grid ( T_tri , t_end_tri , fine_mult =3.0 , n_fine =600 ,
|
|
n_tail =2400)
|
|
# Solve ODE with triangular input
|
|
sol_tri = solve_ivp ( f_rhs , ( t_tri [0] , t_tri [ -1]) , y0 =[0.0 , 0.0] ,
|
|
t_eval = t_tri , args =( vtri_scalar ,) ,
|
|
rtol =1e-8 , atol =1e-11 , method = 'RK45')
|
|
v2_tri = sol_tri.y [1]
|
|
# Numeric area of the triangle ( for normalization )
|
|
vtri_vec = np.array ([ vtri_scalar (t) for t in t_tri ])
|
|
S_tri = np.trapz ( vtri_vec , t_tri )
|
|
# Triangular input & overlay vs analytic impulse response
|
|
plt.figure ( figsize =(8 , 6) )
|
|
plt.subplot (2 ,1 ,1)
|
|
# zoom input on first few T to show the shape
|
|
t_show = np.linspace (0.0 , 5.0* T_tri , 1001)
|
|
v_show = np.array ([ vtri_scalar (t) for t in t_show ])
|
|
plt.plot ( t_show , v_show , 'k')
|
|
plt.grid (True) ; plt.xlabel ( 't (s) ') ; plt.ylabel ( 'V')
|
|
plt.title ( ' Section 2 - Triangular pulse input $v_ { in }(t) $ ')
|
|
|
|
plt.subplot (2 ,1 ,2)
|
|
plt.plot ( t_tri , v2_tri / S_tri , label = r'$v_ {\ mathrm { out }}/\ widehat { S } $(ODE)')
|
|
plt.plot ( t_tri , h (t_tri) , '--' , label = r'$h (t) $ (analytic)')
|
|
plt.grid (True) ; plt.xlabel ( 't (s) ') ; plt.ylabel ( 'V / ( V s ) ')
|
|
plt.title ( ' Triangle : Area - normalized ODE response vs analytic impulse response ')
|
|
plt.legend ( loc = 'best')
|
|
plt.tight_layout () ; plt.show ()
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
# Step input
|
|
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
A_step = 1.0
|
|
t_end_step = 8.0* tau_max
|
|
vstep = lambda t : A_step if t >= 0.0 else 0.0
|
|
# Coarser uniform grid is fine for a step
|
|
t_step = np.linspace (0.0 , t_end_step , 2400)
|
|
# Solve ODE with step input
|
|
sol_step = solve_ivp ( f_rhs , ( t_step [0] , t_step [ -1]) , y0 =[0.0 , 0.0] ,
|
|
t_eval = t_step , args =( vstep ,) ,
|
|
rtol =1e-8 , atol =1e-11 , method = 'RK45')
|
|
v2_step = sol_step.y [1]
|
|
# Analytic step response A * s (t)
|
|
v_step_analytic = A_step * s (t_step)
|
|
# --- PLOT : Step response overlay
|
|
plt.figure ( figsize =(8 , 4.5) )
|
|
plt.plot ( t_step , v2_step , label = ' Step response (ODE)')
|
|
plt.plot ( t_step , v_step_analytic , '--' , label = r' Analytic $A \ , s (t) $ ')
|
|
plt.grid (True) ; plt.xlabel ( 't (s) ') ; plt.ylabel ( 'V')
|
|
plt.title ( ' Step : ODE response vs analytic step response ')
|
|
plt.legend ( loc = 'best')
|
|
plt.tight_layout () ; plt.show ()
|