906-RC回路.sagews

desolve_in_japanese system:sage

Hiroshi TAKEOTO ([email protected])

エンジニアにとって身近な微分方程式をsageを使って解いてみます。

RC回路

RC直列回路の放電

以下のような抵抗（R）とコンデンサー（C）の回路で、

• 抵抗の両端の電圧を$V_R$
• コンデンサーの両端の電圧を$V_C$

とし、$t=0$で$V_C$は$V_0=1$で、すぐに0にした場合のVc電圧の変化$V(t)$を考えてみます。

キルヒホッフの法則から回路を一周したときの電圧降下は０となります。 $$\begin{array}{l l l} V_C + V_R = 0, V_C = V(t) & \cdots & (0) \end{array}$$ となり、コンデンサーからの電流が$i(t)=C\frac{dV(t)}{dt}$、$V_R(t)=Ri(t)$であることから、

$$\begin{array}{l l l} V(t) + RC\frac{dV(t)}{dt} = 0 & \cdots & (1) \end{array}$$ となります。

これをsageを使って解くと以下のようになります。微分方程式の解法にはdesolve関数を使います。

desolve(微分方程式, [求める関数と変数のリスト], [初期値のリスト])

# 使用する変数t, C, Rと関数Vを定義します
var('t C R')
V = function('V',t)
# 微分方程式を定義
de = V + C*R*diff(V,t) == 0
# t=0のV(t)を1として、微分方程式を解く
des = desolve(de,[V,t],[0,1]);view(des)

\newcommand{\Bold}[1]{\mathbf{#1}}e^{-\frac{t}{C R}}

# 解をR=1, C=1のVの変化をプロットする
f(t,R,C) = des
plot(f(t,1,1),[t,0,5])

RC直列回路の充電

つぎに、t=0から電圧V0を掛けて充電した場合のVcの変化を見てみましょう。 以下のような回路で、スイッチを1から2に入れたときの変化になります。

$$\begin{array}{l l l} V(t) + RC\frac{dV(t)}{dt} = V_0 u(t) & \cdots & (2) \end{array}$$ ここで、$u(t)$は、単位ステップ関数、 $$u(t)=\left\{ \begin{array}{l l} 0, & 0 \lt 0 \\ 1, & 1 \ge 0 \\ \end{array} \right.$$ です。

式(1)では右辺が０であり、このような方程式は斉次方程式と呼ばれ、式(2)のように右辺が０でない方程式は非斉次方程式と呼ばれます。 非斉次方程式の解は、

非斉次方程式の一般解　＝　斉次方程式の一般解　＋　非斉次方程式の一つの解（特解）

から求めることができます。

• 非斉次方程式の一つの解（特解）は、スイッチを入れて長い間放っておいた状態を表す解（定常状態の解）であり、特解は$V(t)=V_0$となります。
• desolveの結果から、斉次方程式の一般解は、$V(t) = A e^{\frac{-t}{RC}}$と求まりました。

$t=0$で$V(t)=0$であるとすると、$A=-V_0$となることが分かります。 従って、微分方程式(2)の求める解は、 $$\begin{array}{l l l} V(t) = V_0(1-e^{\frac{-t}{RC}}) & \cdots & (3) \end{array}$$ となります。

v(t, R, C) = (1 - des)
plot(v(t, 1, 1),[t,0,5])

ラブラス変換を使った微分方程式の解法

ラプラス変換を使用して微分方程式の解を求めることができます。ただし、求めることができるのは一般解ではなく「初期値問題」における特殊解です。

ラブラス変換とは

ラプラス変換とは、関数$f(t)$に$e^{-st}$ を掛け、t について0 から無限大まで積分したものです。その積分はs の関数$F(s)$になるが、これをもとの関数$f(t)$のラプラス変換とよび、$\mathcal{L}(f(t))$と表します。 $$F(s) = \mathcal{L}(f(t)) = \int_{0}^{\infty}e^{-st}f(t)dt$$ ここで、$s$は$t$に無関係な変数であり、$t$は実変数である。逆に$f(t)$は$F(s)$の逆変換とよび、$\mathcal{L}^{-1}(F(s))$と表します。

var('t s')
f = function('f', t)
# ラプラス変換のsageでの表現
view(laplace(f, t, s))

\newcommand{\Bold}[1]{\mathbf{#1}}\mathcal{L}\left(f\left(t\right), t, s\right)

一般的なラプラス変換の例

主な関数のラプラス変換をした結果を以下に示します。

元の関数	ラプラス変換結果
$\delta$	1
$1$	$\frac{1}{s}$
$t$	$\frac{1}{s^{2}}$
$t^2$	$\frac{1}{s^{3}}$
$t^{n}$	$s^{{(-n - 1)}} \Gamma\left(n + 1\right)$
$\cos\left(\omega t\right)$	$\frac{s}{{(\omega^{2}+s^{2})}}$
$\sin\left(\omega t\right)$	$\frac{\omega}{{(\omega^{2} + s^{2})}}$
$e^{a t}$	$-\frac{1}{{(a - s)}}$
$t e^{a t}$	$\frac{1}{{(a - s)}^{2}}$

# 変数の宣言と制約条件n > 0をセット
var('n omega a')
assume(n>0)
# 主な関数のラプラス変換の結果
print latex(dirac_delta), laplace(dirac_delta(t), t, s)
print 1, laplace(1, t, s)
print t, laplace(t, t, s)
print t^2, laplace(t^2, t, s)
print latex(t^n), laplace(t^n, t, s)
print latex(cos(omega*t)), laplace(cos(omega*t), t, s)
print latex(sin(omega*t)), laplace(sin(omega*t), t, s)
print latex(exp(a*t)), laplace(exp(a*t), t, s)
print latex(t*exp(a*t)), laplace(t*exp(a*t), t, s)

\delta 1 1 1/s t s^(-2) t^2 2/s^3 t^{n} s^(-n - 1)*gamma(n + 1) \cos\left(\omega t\right) s/(omega^2 + s^2) \sin\left(\omega t\right) omega/(omega^2 + s^2) e^{a t} -1/(a - s) t e^{a t} (a - s)^(-2)

ラプラス変換の性質

ラプラス変換の性質で、特質すべきはラプラス変換の微分によって$s$が掛けられる点です。 $$\mathcal{L}(f') = s\mathcal{L}(f) + f(0)$$ これをsageを使って表現すると以下のようになります。

f = function('f', t)
laplace(diff(f,t), t, s)

s*laplace(f(t), t, s) - f(0)

微分方程式の解法

ラブラス変換を使った微分方程式の解法は、以下の手順で行います。

1. 微分方程式にラプラス変換を施し、微分方程式を変数t領域から変数s領域へ移します
2. 得られた方程式は演算子$s$の代数方程式となるので、代数計算により解を求める
3. 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換する

それでは、sageで微分方程式(1)をラプラス変換を使って解くと、以下のようになります。

# deをラプラス変換する
l1 = laplace(de, t, s); l1

(s*laplace(V(t), t, s) - V(0))*C*R + laplace(V(t), t, s) == 0

# この方程式をlapace(V(t), t, s)について解くと
s1 = solve(l1, laplace(V(t), t, s)); show(s1)
# 解の右辺に初期値V(0)=1を代入する
s11 = s1[0].rhs().subs_expr(V(0) == 1); s11

\left[\mathcal{L}\left(V\left(t\right), t, s\right) = \frac{C R V\left(0\right)}{{(C R s + 1)}}\right]

C*R/(C*R*s + 1)

# 得られた解C*R/(C*R*s + 1)を逆ラプラス変換する
inverse_laplace(s11, s, t)

e^(-t/(C*R))

これで、ラプラス変換からも$e^{\frac{-t}{RC}}$の解を得ることができました。

次に微分方程式(2)に対しても同様にやってみます。

# 微分方程式(2)を定義する
de2 = V + C*R*diff(V,t) == unit_step(t)
l2 = laplace(de2, t, s); l2
# ラプラス変換した後、
s2 = solve(l2, laplace(V(t), t, s)); show(s2)

\left[\mathcal{L}\left(V\left(t\right), t, s\right) = \frac{{(C R s V\left(0\right) + 1)}}{{(C R s^{2} + s)}}\right]

# 得られた解に初期値V(0)=0を代入し
s22 = s2[0].rhs().subs_expr(V(0) == 0); show(s22)
# 逆ラプラス変換する
inverse_laplace(s22, s, t)

\frac{1}{C R s^{2} + s}

-e^(-t/(C*R)) + 1

微分方程式(1)と同様微分方程式(2)も期待した解を得ることができました。

maximaを使った微分方程式の解法

sageでもある程度の微分方程式は解くことができますが、初期値の柔軟性という点では、maximaに一日の長があるようです。

次に以下の微分方程式をmaximaを使って解いてみます。 $$x^2\frac{dy}{dx}+3xy = \frac{sin(x)}{x}, y(\pi) = 0$$ maximaを使った微分方程式の解法は以下の手順で行います。

• maximaインタフェースを使って微分方程式を定義する
• ode2関数を使って微分方程式を解く
• ic1またはic2を使って初期条件をセットする

_sage_()関数は、別の処理系（maxima, scilab, R等）の結果をsageの表現に変換するために使用します。

# 変数x, yを定義
var('x y')
# maximaを使って微分方程式を定義：diffの前の'に注意
deq = maxima("x^2*'diff(y,x) + 3*x*y = sin(x)/x")

# 微分方程式を解く
d2 = deq.ode2(y,x).ic1(x=pi,y=0)
# 解を表示
show(deq.ode2(y,x))
# maxima の結果をsageの表現に変換するには_sage()_を使い、rhs()を使ってyの関数のみを取り出すことができる
d3 = d2._sage_(); d3.rhs()

y={{{\it c}-\cos x}\over{x^3}}

-(cos(x) + 1)/x^3

# ic1を使って初期値を指定して微分方程式を解く
d4 = deq.ode2(y,x).ic1(x=pi,y=0); d4

y=-(cos(x)+1)/x^3

# 結果をプロットする
plot(d4._sage_().rhs(), [x, 0, pi])

このようにsageを使って微分方程式を様々な手法で解くことにより、解の検証も可能です。

是非、この機会にsageを使って微分方程式を解いてみてください。