ラプラス変換をやろう
RLC回路
交流RLC回路
１。ラプラス変換のやり方
このワークシートはMath by Codeの一部です。
前回まで、常微分方程式のタイプ分けをタイプに応じた解法を見てきた。
今回からは、視点を変えて解法の種類を広げてみよう。まずはラプラス変換。
＜ラプラス変換＞
ラプラス変換は、微分方程式を代数式に変えて、代数的に単純化する、
それを逆変換するだけで、微分方程式が解けてしまう。
つまり、微分積分計算を代数計算に置き換えて、手計算を減らすことができるというおいしい手順の１つだ。ラプラスさんの編み出した方法は、電気技師のヘビサイドという人物が演算子を使って電気回路の設計に使われて広がったといわれている。制御システムを作るのにも使われる。
もちろんどんな微分方程式でも解けるというわけではないにしろ幅広いタイプに使える。
・もとの関数f(t)に対するラプラス変換 L(f)=integral(exp(-st) f(t), dt, [0, ∞])= ∫e-st f(t) dt
不定積分ではなく定積分の形になったいるが、[0, ∞]が固定なので、最後の数値代入以外は省略する。
たて棒の上下に積分範囲をかくスタイルもあるが、行の節約から[∞,0]閉区間のように記すことにする。
積分範囲に∞があるので特異な積分だけど、数Tを入れてからTを∞にすると考えればよいね。
積分変数はtだが、指数関数の変数がtの-s倍になっていることに着目しよう。
指数関数exp(-st)をf(t)にかけることで、f(t)を急速に0に近づけることで、片側無限区間の積分値が
収束するようにしているね。
・ラプラス変換は積分の１種なので線形性がある。だから、基本となる関数の変換公式の組み合わせで
変換ができてしまう。逆変換も同様。
(例)
L(1)=∫e-st 1dt= -1/s e-st[∞,0]＝0-(-1/s)=1/s＝F(s)
L(eat)=∫e(a-s)tdt= 1/(a-s) e(a-s)t[∞,0]＝0-(1/(a-s))=1/(s-a)＝F(s)
ステップ関数θ(非負で１、負で０）L(θ)=∫e-st θdt=∫e-st 1dt=1/s
デルタ関数(δ(t)=0 , ∫g(t)δ(t)[∞,-∞] dt = g(0) ) L(δ)=∫e-st δ(t)dt=e-s 0 =1
・第一推移定理L(f(t))=F(s)ならば、L(eatf(t))=F(s-a) ,L(e-atf(t))=F(s+a)
・相似定理 L(f(t))=F(s)ならば、aが正のときL(f(at))=1/aF(s/a)
・基本変換集( f(t)→L(f))
べき関数1→1/s , t → 1/s2, t2 →2!/s3, nが整数のとき 、aが正数のとき
指数関数eat→ (aが複素数でも）
三角関数cosωt→ , sinωt→ ,cosh at=(eat+e-at)/2→, sinh at=(eat-e-at)/2→
振幅変化三角関数 eat cosωt→ , eat sinωt→
・微分の変換集 (ｓはあるkより大きい、そこではfは連続関数）初期値f(0)=zとする。f→Yとする。
微分を1階するごとにYにsを1回かける。
L(f')=sY-z （sの1次式)
L(f'')=s2Y-sz-z' （sの2次式）
....................................... （n解のLはｎ次式だけど、最高次以外はzのn-1階以下の微分の初期値が係数）
L(f(n))=snY-(sn-1z+sn-2z'+......+ z(n-2)) sの指数とｚの微分階数の和がn-1のsのn-1次多項式をひく。
・積分の変換(微分のs倍の逆だからsで割る）
L(∫f(r) dr, [0,t] ) = 1/s Y
・部分分数展開の例
・ラプラス変換の手順
（ステップL）
微分方程式y''+ay'+by=r(t) , 初期値z=K0, z'=K1とする。L(y)=Y, L(r)=Rとすると、
1y''+ay'+by→1 L( y'')+aL(y')+bL(y)=1(s2Y-sz-z' )+a(sY-z )+b Y=(1 s2+as+b) Y - (s+a)z- z' =R
（ステップY）
Q= 1/ (s2+as+b)とおくと、Y={(s+a)z+z'} Q+RQ。z=z'=0なら、Y=RQとなりQ=Y/R＝L(出力)/L(入力）
Yを部分分数展開などを使って、基本変換集のL(f)にあたる式の和に分解する。
（ステップL-1)
基本変換集のL(f)部分をf(t)に逆変換する。
（例）
y''-y=t , f(0)=zとするとき、z=1, z'=1。L(y)=Yとする。
（Lする) Y(s2-1)=sz+z'+1/s2= s+1+1/s2
（Yでとく）Q=1/(s2-1)とおく、Y=(s+1+1/s2)Q=(s+1)/(s2-1)+1/s2(s2-1)=1/(s-1)+(1/(s2-1)-1/s2)
（逆Lする）f(t)=L-1(Y)=L-1[1/(s-1)]+L-1[1/(s2-1)]-L-1[1/s2]=et+sinh t -t =et+(et-e-t)/2-t=3/2et-1/2e-t -t
2。ラプラス変換になれよう。
（例）
yはt秒後の自由落下の高さ
加速度はy'’=-g, 最初の高さf(0)=y0, 最初の速度f'(0)=0のとき
y''=-g , f(0)=zとするとき、z=y0, z'=0。L(y)=Yとする。
（Lする) Ys2-(sz+z')=-g/s Ys2-sz=-g/s Ys2=sz-g/s
（Yでとく）Q=1/s2とおく、Y=(sz-g/s)Q=z/s-g/s3
（逆Lする）f(t)=L-1(Y)=zL-1[1/s]-gL-1[1/s3]=z - 1/2g t2=y0- 1/2g t2
（例）
y'+2y=e-t, z=z'=0
（Lする) Y(s+2)=z+1/(s+1)=1/(s+1)
（Yでとく）Q=1/(s+2)おく、Y=1/(s+1)(s+2)=(1/(s+1)- 1/(s+2))
（逆Lする）f(t)=L-1(Y)=zL-1[1/(s+1)]-L-1[1/(s+2)]t=e -t- e-2t
（例)
y''+3y'+2y=eit, z=z'=0
（Lする) Y(s2+3s+2)=1/(s-i)
特性方程式t2+ 3s+2=0の解はp,q=-1,-2だから、Y(s+1)(s+2)=1/(s-i)
（Yでとく） Y=1/(s-i)(s+1)(s+2)=
=
=
（逆Lする）y(t)=
=
（例）
y''+3y'+2y=sin t , z=z'=0 （A)
y''+3y'+2y=cos t , z=z'=0（B)
Aは上の解の虚部Im(y(t))、Bは上の解の実部Re(y(t))
３．実装
質問：ラプラス変換をコードやるにはどうしたらよいでしょうか。
＜Pytho with Sympy＞
Sympyをつかいます。
残念ですがワンストップでは、計算できません。
1回目の入力で、ラプラス変換にはtの関数式f(t)を入力して渡します。
ラプラス変換の結果はsympy流の記述で等式が出力されます。
2回目の入力で、結果の等式を目で見ながら、通常に近い形の等式を入力します。
Yで解き、部分分数計算し、逆変換するを一気にできます。
青い部分が毎回変わる入力と出力です。他は、入力内容に応じて変わることがあるかもしれません。
Eq(左辺,右辺)で等式を作ります。laplace_correspondenceで対応する関数{f: F}を等式に伝えましょう。
F = F.subs({z:1, dz:1})は等式Fに初期値を代入します。
solve(F,Y)は等式FをY＝にした右辺を答えます。
apart(s)はsの有理関数を部分分数に分解します。
手動のプロセスと比べてみよう。
y''-y=t , f(0)=zとするとき、z=1, z'=1。L(y)=Yとする。
（Lする) Y(s2-1)=sz+z'+1/s2= s+1+1/s2
（Yでとく）Q=1/(s2-1)とおく、Y=(s+1+1/s2)Q=(s+1)/(s2-1)+1/s2(s2-1)=1/(s-1)+(1/(s2-1)-1/s2)
（逆Lする）f(t)=L-1(Y)=L-1[1/(s-1)]+L-1[1/(s2-1)]-L-1[1/s2]=et+sinh t -t =et+(et-e-t)/2-t=3/2et-1/2e-t -t
#{IN]part1
from sympy import *
from sympy.abc import t, s
y = Function('y')(t) #もとの関数
Y = Function('Y')(s) #ラプラス変換後の関数
ddy=Derivative(y, t, 2)
dy=Derivative(y, t)
f_left = ddy - y
f_right = t
ode = Eq(f_left, f_right)
print(f"もと: {ode}")
F_left = laplace_transform(f_left, t, s, noconds=True)
F_right = laplace_transform(f_right, t, s, noconds=True)
lap = Eq(F_left, F_right)
F = laplace_correspondence(lap, {y : Y})
print(f"Lのあと:{F}")
#[OUT]もと: Eq(-y(t) + Derivative(y(t), (t, 2)), t)
Lのあと:Eq(s**2*LaplaceTransform(y(t), t, s) - s*y(0) - LaplaceTransform(y(t), t, s) - Subs(Derivative(y(t), t), t, 0), s**(-2))
#[IN]part2
#手動で、LaplaceTransform(y(t), t, s)をY、y(0)をｚ、Subs(Derivative(y(t), t), t, 0)をdzとして
#ラプラス変換後の式をかく。初期値z=y(0)=1,dz=y'(0)=1を後でsubstituteする。
z, dz = symbols('z, dz')
F = Eq(s**2*Y - s*z - Y - dz, s**(-2))
F = F.subs({z:1, dz:1})
print(f"変換後の等式:{F}")
#Yで解いて、部分分数分解して、逆L
Y_sol = solve(F, Y)[0]
Y_sol_apart = Y_sol.apart(s)
print(f"Y = {Y_sol_apart}")
y_t = inverse_laplace_transform(Y_sol_apart, s, t)
y_t = y_t.subs({Heaviside(t):1})
print(f"y ={y_t}")
#[OUT]
変換後の等式:Eq(s**2*Y(s) - s - Y(s) - 1, s**(-2))
Y = -1/(2*(s + 1)) + 3/(2*(s - 1)) - 1/s**2
y =-t + 3*exp(t)/2 - exp(-t)/2
質問：複素関数のラプラス変換をコードやるにはどうしたらよいでしょうか。
＜Pytho with Sympy＞
pythonでは虚数単位がiではなくjになることに気をつければ、実数と同様に作動します。
手動のプロセスと比べてみよう。
#y''+3y'+2y=e^it, z=z'=0
from sympy import *
from sympy.abc import t, s
y = Function('y')(t) #もとの関数
Y = Function('Y')(s) #ラプラス変換後の関数
ddy=Derivative(y, t, 2)
dy=Derivative(y, t)
f_left = ddy + 3*dy + 2*y
f_right = exp(j*t)
ode = Eq(f_left, f_right)
print(f"もと: {ode}")
F_left = laplace_transform(f_left, t, s, noconds=True)
F_right = laplace_transform(f_right, t, s, noconds=True)
lap = Eq(F_left, F_right)
F = laplace_correspondence(lap, {y : Y})
print(f"Lのあと:{F}")
#[OUT]
もと: Eq(2*y(t) + 3*Derivative(y(t), t) + Derivative(y(t), (t, 2)), exp(j*t))
Lのあと:Eq(s**2*LaplaceTransform(y(t), t, s) + 3*s*LaplaceTransform(y(t), t, s) - s*y(0) + 2*LaplaceTransform(y(t), t, s) - 3*y(0) - Subs(Derivative(y(t), t), t, 0), 1/(-j + s))
#（Lする) Y(s2+3s+2)=1/(s-i)#手動で、LaplaceTransform(y(t), t, s)をY、y(0)をｚ、Subs(Derivative(y(t), t), t, 0)をdzとして
#ラプラス変換後の式をかく。初期値z=y(0)=0,dz=y'(0)=0を後でsubstituteする。
z, dz = symbols('z, dz')
F = Eq(s**2*Y + 3*s*Y -s*z + 2*Y -3*z -dz, 1/(-j +s))
F = F.subs({z:0, dz:0})
print(f"変換後の等式:{F}")
#Yで解いて分数分解して逆L
Y_sol = solve(F, Y)[0]
Y_sol_apart = Y_sol.apart(s)
print(f"Y = {Y_sol_apart}")
y_t = inverse_laplace_transform(Y_sol_apart, s, t)
y_t = y_t.subs({Heaviside(t):1})
print(f"y ={y_t}")
#[OUT]
変換後の等式:Eq(s**2*Y(s) + 3*s*Y(s) + 2*Y(s), 1/(-j + s))
Y = 1/((j + 2)*(s + 2)) - 1/((j + 1)*(s + 1)) + 1/((-j + s)*(j + 1)*(j + 2))
y =exp(-2*t)/(j + 2) - exp(-t)/(j + 1) + exp(j*t)/((j + 1)*(j + 2))
複素数の計算を最低限しかやってくれないようなので、違う答えのように見えてしまいますが同じです。
y(t)=
＜geogebra＞
geogebraでもラプラス変換、逆変換の関数はあるようですが、
基本関数が計算できるだけで、微分の変換ができないようです。だから、微分方程式を解く用途には
使えないですね。
４．実例
コンデンサーは文字通り、電気エネルギー（電荷）を蓄える。
さらに、直流を通さず、交流は周波数が高いほど通しやすい。
抵抗Rだけでなく、コンデンサーCをつないだ回路で、電流変化を考えてみます。
＜RC回路＞
RC回路はRとCに電圧V0をかけた回路です。
キルヒホフ第2法則（電池で高まった電圧はCとRで下がり、電池のマイナス側で０に戻る）を使います。
初期条件Q(0)=0でC,R,V0を順につないでスイッチを入れたときの回路に流れる電荷Qを求めよう。
Qは時間tの関数とします。電流Iは電荷Qの時間変化なので、I=dQ/dtです。
コンデンサーの両端の電圧V=Q/Cで、抵抗Rの両端の電圧V=IR＝RQ'
電圧の和がV0に等しいから、R Q' + Q/C = V0 式変形 Q'+Q/CR=V0/R
見やすさのために、一時的にQをyとして表記し、最後にQに戻すことにします。
y'+y/CR=V0/R , y(0)=0とするとき、z=0。L(y)=Yとする。
（Lする) Y(s+1/CR)=z+V0/R/s=V0/R/s
（Yでとく）Y=(V0/R){1/s*1/(s+1/CR)}=(V0/R){1/s-1/(s+1/CR)}CR=CV0{1/s-1/(s+1/CR)}
（逆Lする）f(t)=L-1(Y)=CV0{L-1[1/s]-L-1[1/(s+1/CR]}=CV0(1 - e-t/CR)
Q=CV0(1 - e-t/CR)
e-t/CRはt=0あたりではほぼ１でｔが∞で０に近づく。Qは０から増えてCV0で飽和する。
また、I=dQ/dt=(CV0(1 - e-t/CR))'=CV0 (-e-t/CR)'=CV0/CR e-t/CR=(V0/R) e-t/CR
e-t/CRはt=0あたりではほぼ１でｔが∞で０に近づく。Iは(V0/R)から減って0になる。
RC回路
＜直流RLC回路＞
Lはコイルです。RLC回路は、RC回路にさらにコイルをつないだものです。
回路といっても、計算上のものなので、数式の上での前提での架空の推論をします。
工作をして実験するわけではありません。
それぞれの電圧は、V=IR,V=L dI/dt, V=Q/Cとします。
定数はR=2、L=1, C=0.5, V0=10とし、
コンデンサーの電荷Q, 回路の電流Iは時間tの関数で、初期値Q(0)=I(0)=0としましょう。
L(I)=Yとします。電流Iは電荷Qの時間変化なので、I=dQ/dt=Q'です。 L(Q)= 1/s L(I)=1/s Y。
z=z'=0
この回路に減衰しながら振動する電流が流れるのはどんなときでしょう。
R,L,Qの電圧の和はRI+L dI/dt + Q/C=V0
（Lする）RY+sLY+1/CsY=V0/s sR/LY+s2Y+1/CL Y=V0/L
（Yでとく） Y=V0/L(s2+R/L s +1/CL)
特性方程式 t2+R/L s + 1/CL=0の判別式D=(R/L)2-4/CL=(1/L)2(R2-4L/C)が負のとき、
（つまり、Lに対してC,Rが小さいとき、C容量・抵抗RがインダクタンスLに対して小さいとき）
(t- α)2+β2=0と平方完成されるとして , V0/L=pとおくと、
Y=p *1/{(s- α)2+β2}=p/β *β/{(s- α)2+β2}
（逆Lする）I=p/β *eαt sinβt α=-R/2L , β=√-(D/4)
定数を入れると、p=10/1=10, R/L=2/1=2, 1/CL=1/(1*0.5)=2,
α=-R/2L=-2/2*1=-1,D=2^2-4*2=-4, β=√-(-4/4)=1
I(t) =10/ 1 * e-t sin t = 10 e-t sin t
直流RLC回路の減衰電流
＜交流RLC回路＞
回路の電圧が時間とともに変化できる、つまり交流電圧Vsinωtをかけた回路ではどうでしょうか。
RI+L dI/dt + Q/C=Vsinωtです。
複素数の範囲で計算することで、複素数の範囲で解きましょう。
R,L,Qの電圧の和はRI+L dI/dt + Q/C=Veiωtとします。（iω=c, V/L=Aとおく）c2=-ω2
L(Q)=Yとします。I=dQ/dt=Q'です。z=z'=0
RI+L dI/dt + Q/C=Vect
L Q'' +R Q'+ 1/C Q=Vect
（Lする）Y(s2+R/L + 1/CL)=A/(s-c)
（Yでとく） Y=A/(s2+R/L s +1/CL)(s-c)
特性方程式 t2+R/L s + 1/CL=0の2解をa,bとすると、(t-a)(t-b)=0
a=-R/2L ＋i√D/2L , b = -R/2L -i√D/2L a,bの実部は負の数
a+b=-R/L, ab=1/CL
Y=A/(s-a)(s-b)(s-c)=A*/(a-b){1/(s-a)(s-c)-1/(s-b)(s-c)}
=A/(a-b)[1/(a-c) * {1/(s-a)-1/(s-c)}- 1/(b-c)*{1/(s-b)-1/(s-c)}]
= A/(a-b)[1/(a-c)(s-a)-1/(a-c)(s-c)-1/(b-c)(s-b)+1/(b-c)(s-c)]
=-A/(a-b)[1/(c-a)(s-a)+1/(b-c)(s-b)+1/(s-c){1/(a-c)-1/(b-c}]
=-A/(a-b)[1/(c-a)(s-a)+1/(b-c)(s-b)+1/(s-c){(b-a)/(a-c)(b-c)]
=-A/(a-b)[(b-c)/(b-c)(c-a)(s-a)+(c-a)/(c-a)(b-c)(s-b)+(a-b)/(s-c)(c-a)(b-c)]
=-A/{(a-b)(b-c)(c-a)[(b-c)/(s-a)+(c-a)/(s-b)+(a-b)/(s-c)]
=
（逆Lする）Q(t)=
a,bの実部は負の数だからtを十分大きくすると、eat, ebtの極限は０になる。
Q(t)の極限は
=(V/L)/[-ω2+R/L iω+1/CL] eiωt
= V/[-Lω2+R iω+1/C] eiωt
ラプラス変換で求めた解から減衰する項を消し、交流定常状態の解を出すことができました。
この定常状態の虚部を求めてみよう。
分母Z(ω)=−Lω2+Riω+1/C1=(1/C−Lω2)+i(Rω)=A+iB (複素インピーダンスという)と置くと、
|Z(ω)|=A2+B2で、その位相ϕ=arg(Z(ω))=arctan(B/A)です。
だから、Q(t)=Veiωt/(∣Z(ω)∣eiϕ)=Vei(ωt−ϕ)/|Z(ω)|=V/|Z(ω)| * (cos(ωt−ϕ)+i sin(ωt−ϕ))
この Q(t) の虚部は、Im[Q(t)]=V sin(ωt−ϕ)/|Z(ω)|