PythonでRC回路の過渡解析をしてみる その2 数値計算 by Scipy編

前回建てた微分方程式と初期条件
Vc'=1/(RC)(Vin-Vc)

Vin(0) = 0.0
Vc(0)= Vcc

前回、手計算で解いた結果
Vc(t)=Vcc×exp{-1/(RC)×t}

でも、定数変化法で解くので手間が掛かる。


手計算ではなく数値計算する方法。

Pythonの数値計算モジュールScipyを使う。図１に実行結果を示す。

図１：実行結果

from scipy import *

from pylab import *

deriv = lambda y,t : array([y[1],-y[0]])

# Initial parameters

R=10**3

C=10**-7

Vcc=5.





def Vin(t):

    return 0.



def Vcdash(Vc,t):

    return (Vin(t)-Vc[0])/(R*C)



start=0

end=10**-3

numsteps=10000

time=linspace(start,end,numsteps)



from scipy import integrate

Vc0=Vcc

Vc=integrate.odeint(Vcdash,Vc0,time)



#Plot Setting

figure(1,figsize=(6,4)).subplots_adjust(bottom=0.11)

xlabel('t')

grid(True)

hold(True)



plot(time,Vc)

#Setting Window title and Save Figure

title('Transient Analysis of RC Circuit')

savefig('transient analysis.png',dpi=72)



show()

integrate.odeint()関数には、微分方程式と初期値と時間に相当する行列を渡すと計算してくれます。
微分方程式はこの場合Vcdash()です。
なんで微分方程式が
def Vcdash(Vc,t):
return (Vin(t)-Vc[0])/(R*C)
って表現されるのか？その理由を以下に示します。

＜リストを用いて表現される理由＞
微分方程式を数値解析的に解くには差分方程式だったか漸化式だったかを用いる。（たしか）
差分方程式@wikipedia

解きたい微分方程式

Vc'=1/(RC)(Vin-Vc)

を前進差分で差分方程式にしてやると、
Vc'=(Vc(t+h)-Vc(t))/h=1/(RC)(Vin-Vc(t)) --(1)
ここで、hを数値計算するメッシュの分解能として考えると（＝メッシュがhずつ刻まれていると考えると）、
Vc(t) =V[n]
Vc(t+h)=v[n+1]
とおける。
すると(1)式は

(Vc[t+1]+V[t])/h=(Vin-Vc[t])/RC
V[t+1]=Vin/RC-Vc[t](h-1/RC)

となり、ここでhを0に近づければ、（計算機では０にはできないが）
Vc[t+1]≒1/(RC)(Vin-Vc[t]) --(2)


(1),(2)を見比べると、Vc'はVc[t+1]に、VcはVc[t]に置き換わっただけということが分かる。ということで、微分方程式はリストを使って表現できる。
2階微分方程式の場合、例えば、

y''=ay'+by+c

を解きたい場合には、

y'=x

とおくと

x'=ax+by+c

のように、1階常微分方程式になりscipy.integrate.odeint()関数で扱える形式になる。
ただし、今度は、yの初期値だけでなくy'=xの初期値も必要となってくるので注意が必要。
参考Scipy CookBook




ついでに、このプログラムのVin(t)という関数は電源電圧を表現している、これをsin(ωt)を返す関数にした例↓。



図２:Vin(t)=Vcc*sin(10**5*t)とした場合の実行結果