Skip to content

Instantly share code, notes, and snippets.

/gist:1190103
Created Sep 2, 2011

Embed
What would you like to do?
################################################################################
# undriven RLC circuit
# solved for analytical solution in Maxima
# see these MIT OCW notes for a good intro:
# http://ocw.mit.edu/courses/physics/8-02sc-physics-ii-electricity-and-magnetism-fall-2010/undriven-rlc-circuits/
# Author: Josh Stults
# Date: 28 Aug 2011
# Website: www.variousconsequences.com
#
# governing equations as defined in Maxima:
# depends(I,t);
# gov_eqn : expand((diff(I,t,2) * L + R*diff(I,t) + I/C = 0)/L);
# sol_pos : ode2(gov_eqn, I, t); /* C(CR^2-4L) > 0 */
# sol_neg : ode2(gov_eqn, I, t); /* C(CR^2-4L) < 0 */
# sol_zero : ode2(gov_eqn, I, t); /* C(CR^2-4L) == 0 */
# sol_pos_ic : ic2(sol_pos, t=0, I=0, diff(I,t)=v0/L);
# sol_neg_ic : ic2(sol_neg, t=0, I=0, diff(I,t)=v0/L);
# sol_zero_ic : ic2(sol_zero, t=0, I=0, diff(I,t)=v0/L);
import scipy as sp
from scipy.integrate import ode
from scipy.optimize import fmin
bootindex = sp.random.random_integers
import matplotlib.pylab as plt
# dump these expressions from Maxima using f90; then add 'sp.' and remove '&'
def undriven_rlc(t,v0,R,L,C):
check = C*(C*R**2 - 4*L)
if(check>0):
I = v0*sp.exp(t*(sp.sqrt(R**2/L**2-4/(C*L))-R/L)/2.0)/(L*sp.sqrt(R**2/L**2-4/(C*L)))-v0*sp.exp(t*(-sp.sqrt(R**2/L**2-4/(C*L))-R/L)/2.0)/(L*sp.sqrt(R**2/L**2-4/(C*L)))
elif(check==0):
I = t*v0*sp.exp(-t*R/L/2.0)/L
else:
I = 2*v0*sp.exp(-t*R/L/2.0)*sp.sin(t*sp.sqrt(4/(C*L)-R**2/L**2)/2.0)/(L*sp.sqrt(4/(C*L)-R**2/L**2))
return(I)
def undriven_rlc_ls(X,v0,t1,i1,t2,i2):
R = X[0]
L = X[1]
C = X[2]
err_1 = undriven_rlc(t1,v0,R,L,C) - i1
err_2 = undriven_rlc(t2,v0,R,L,C) - i2
return(sp.sqrt(sp.sum(sp.append(err_1,err_2)**2)))
################################################################################
v0 = 800.0 # V initial potential on the capacitor
R = 39e-3 + 1e-3 + 33e-3 + 227e-3 # Ohms
C = 3e-4 # F
L = (C*R**2)/4.0 #6.75e-6 # H
X0 = sp.array([R,L,C])
print (C*R**2)/4.0
print C*(C*R**2 - 4*L)
t = sp.linspace(0,1e-3,4000)
I0 = undriven_rlc(t,v0,R,L,C)
d = 0.25
Ip = undriven_rlc(t,v0,R,L/d,C)
Im = undriven_rlc(t,v0,R,L*d,C)
# read in g3data scraped current traces
t1,i1,t1err,i1err = sp.loadtxt('IMG_1663.JPG.dat',unpack=True)
t2,i2,t2err,i2err = sp.loadtxt('IMG_1663.JPG_2.dat',unpack=True)
toff = t1.min()
t1 = t1 - toff
t2 = t2 - toff
# fit least squares estimate to the data
xopt = fmin(undriven_rlc_ls, X0,(v0,t1,i1,t2,i2))
tfit = sp.linspace(0,0.0025,4000)
ifit = undriven_rlc(tfit,v0,xopt[0],xopt[1],xopt[2])
# take a look at the autocorelation of the residuals
resid1 = i1 - undriven_rlc(t1,v0,xopt[0],xopt[1],xopt[2])
fresid1 = sp.fft(resid1)
autocor1 = sp.ifft(fresid1 * fresid1.conj()).real
autocor1 = autocor1/autocor1.max()
resid2 = i2 - undriven_rlc(t2,v0,xopt[0],xopt[1],xopt[2])
fresid2 = sp.fft(resid2)
autocor2 = sp.ifft(fresid2 * fresid2.conj()).real
autocor2 = autocor2/autocor2.max()
# bootstrap the residuals to give confidence intervals
nboots = 100
xboot = sp.zeros((nboots,3))
for i in xrange(nboots):
bi1 = bootindex(0,i1.shape[0]-1,i1.shape[0])
bi2 = bootindex(0,i2.shape[0]-1,i2.shape[0])
xboot[i] = fmin(undriven_rlc_ls, xopt, (v0,t1,i1+resid1[bi1],t2,i2+resid2[bi2]))
plt.figure()
plt.plot(t, I0, label='L=%g H'%L)
plt.plot(t, Ip, label='L=%g H'%(L/d))
plt.plot(t, Im, label='L=%g H'%(L*d))
plt.xlabel('t (s)')
plt.ylabel('I (A)')
plt.legend(loc=0)
ax = plt.gca()
plt.text(0.5,0.5,'initial voltage: %g V\nresistance: %g Ohms\ncapacitance: %g F'%(v0,R,C),fontsize='large',transform=ax.transAxes)
plt.savefig("undriven_rlc_pulse.png")
plt.figure()
plt.plot(t1,i1,'o')
plt.plot(t2,i2,'o')
plt.plot(tfit,ifit,'k-')
ax = plt.gca()
plt.text(0.2,0.6,'Given\n initial voltage: %g V\n IMG_1663.JPG scraped w/g3data\nFit\n resistance: %1.2e [%1.2e,%1.2e] Ohms\n inductance: %1.3g [%1.2e,%1.2e] H\n capacitance: %1.2e [%1.2e,%1.2e] F'%(v0,xopt[0],xboot[:,0].min(),xboot[:,0].max(),xopt[1],xboot[:,1].min(),xboot[:,1].max(),xopt[2],xboot[:,2].min(),xboot[:,2].max()),fontsize='large',transform=ax.transAxes)
plt.xlabel('t (s)')
plt.ylabel('I (A)')
plt.legend(loc=0)
plt.savefig("undriven_rlc_fit.png")
plt.figure()
# significant autocorrelation?
plt.plot(autocor1[0:autocor1.shape[0]/2],'o')
plt.plot(autocor2[0:autocor2.shape[0]/2],'o')
plt.ylabel('autocorrelation')
plt.xlabel('lag')
plt.savefig("fit_resid_autocor.png")
plt.figure()
# look about normal?
nbins = 24
plt.hist(resid1,bins=nbins,normed=True,alpha=0.3)
plt.hist(resid2,bins=nbins,normed=True,alpha=0.3)
plt.savefig("fit_resid_hist.png")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.