sage/微分方程式を解いてみる
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2009/11/01からのアクセス回数 &counter;
sageを使ってエンジニアにとって身近な微分方程式を解いてみ...
また、実際にsageを使って出力されたワークシートは、以下のU...
http://www.sagenb.org/home/pub/906/
** RC回路 [#f65f28fa]
*** RC直列回路の放電 [#u93cc6d1]
以下のような抵抗（R）とコンデンサー（C）の回路((図と式の...
&ref(RC1.gif);
- 抵抗の両端の電圧をVr
- コンデンサーの両端の電圧をVc
とし、t=0でVcにv0の電圧を掛け、すぐに0にした場合のVc電圧...
キルヒホッフの法則から回路を一周したときの電圧降下は０と...
$$
\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}\...
$$
\begin{array}{l l l}
V(t) + RC\frac{dV(t)}{dt} = 0 & \cdots & (1)
\end{array}
$$
となります。
これをsageを使って解くと以下のようになります。微分方程式...
desolve(微分方程式, [求める関数と変数のリスト], [初期値...
sageへの入力：
#pre{{
# 使用する変数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)
}}
&color(blue){\(e^{-\frac{t}{C R}}\)};
#pre{{
# 解をR=1, C=1のVの変化をプロットする
f(t,R,C) = des
plot(f(t,1,1),[t,0,5])
}}
&ref(sage0.png);
*** RC直列回路の充電 [#z90b8582]
つぎに、t=0から電圧\(V_0\)を掛けて充電した場合のVcの変化...
&ref(RC2.gif);
$$
\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)では右辺が０であり、このような方程式は斉次方程式と呼...
非斉次方程式の一般解 ＝ 斉次方程式の一般解 ＋ 非斉次...
から求めることができます。
- 非斉次方程式の一つの解（特解）は、スイッチを入れて長い...
- desolveの結果から、斉次方程式の一般解は、\(V(t) = A e^{...
t=0でV(t)=0であるとすると、A=−V0 となることが分かります。...
$$
\begin{array}{l l l}
V(t) = V_0(1-e^{\frac{-t}{RC}}) & \cdots & (3)
\end{array}
$$
となります。
sageへの入力:
#pre{{
v(t, R, C) = (1 - des)
plot(v(t, 1, 1),[t,0,5])
}}
&ref(sage0-1.png);
** ラブラス変換を使った微分方程式の解法 [#aa81fa10]
ラプラス変換を使用して微分方程式の解を求めることができま...
*** ラブラス変換とは [#ta611804]
ラプラス変換とは、関数\(f(t)\)に\(e^{-st}\) を掛け、t に...
$$
F(s) = \mathcal{L}(f(t)) = \int_{0}^{\infty}e^{-st}f(t)dt
$$
ここで、\(s\)は\(t\)に無関係な変数であり、\(t\)は実変数で...
sageへの入力:
#pre{{
var('t s')
f = function('f', t)
# ラプラス変換のsageでの表現
view(laplace(f, t, s))
}}
&color(blue){\(\mathcal{L}\left(f\left(t\right), t, s\rig...
*** 一般的なラプラス変換の例 [#u005628c]
主な関数のラプラス変換をした結果を以下に示します。
|元の関数|ラプラス変換結果|h
|\(\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(\omega t)\)|\(\frac{s}{{(\omega^{2}+s^{2})}}\)|
|\(sin(\omega t)\)|\(\frac{\omega}{{(\omega^{2} + s^{2})}...
|\(e^{a t}\)|\(-\frac{1}{{(a - s)}}\)|
|\(t e^{a t}\)|\(\frac{1}{{(a - s)}^{2}}\)|
上記の表を出力するsageの入力：
#pre{{
# 変数の宣言と制約条件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 t^n, laplace(t^n, t, s)
print cos(omega*t), laplace(cos(omega*t), t, s)
print sin(omega*t), laplace(sin(omega*t), t, s)
print exp(a*t), laplace(exp(a*t), t, s)
print t*exp(a*t), laplace(t*exp(a*t), t, s)
}}
出力：
#pre{{
\delta 1
1 1/s
t s^(-2)
t^2 2/s^3
t^n s^(-n - 1)*gamma(n + 1)
cos(omega*t) s/(omega^2 + s^2)
sin(omega*t) omega/(omega^2 + s^2)
e^(a*t) -1/(a - s)
t*e^(a*t) (a - s)^(-2)
}}
*** ラプラス変換の性質 [#bd46fa5d]
ラプラス変換の性質で、特質すべきはラプラス変換の微分によ...
$$
\mathcal{L}(f') = s\mathcal{L}(f) + f(0)
$$
これをsageを使って表現すると以下のようになります。
sage入力：
#pre{{
f = function('f', t)
laplace(diff(f,t), t, s)
}}
&color(blue){s*laplace(f(t), t, s) - f(0)};
*** 微分方程式の解法 [#pdc59ce1]
ラブラス変換を使った微分方程式の解法は、以下の手順で行い...
+ 微分方程式にラプラス変換を施し、微分方程式を変数t領域か...
+ 得られた方程式は演算子$s$の方程式となり、solve関数によ...
+ 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換...
それでは、sageで微分方程式(1)をラプラス変換を使って解くと...
sageの入力：
#pre{{
# deをラプラス変換する
l1 = laplace(de, t, s); l1
}}
&color(blue){(s*laplace(V(t), t, s) - V(0))*C*R + laplace...
#pre{{
# この方程式を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
}}
&color(blue){(laplace(V(t), t, s) == C*R*V(0)/(C*R*s + 1)};
&color(blue){C*R/(C*R*s + 1)};
#pre{{
# 得られた解C*R/(C*R*s + 1)を逆ラプラス変換する
is1 = inverse_laplace(s11, s, t); is1
}}
&color(blue){e^(-t/(C*R))};
これで、ラプラス変換からも\(e^{\frac{-t}{RC}}\)の解を得る...
手計算では、ラプラス変換表を使って逆ラプラス変換を行って...
- \(\frac{CR}{(CRs + 1}\)を\(-CR\)で割ると\(\frac{-1}{(\f...
- ラプラス変換表の\(e^{at}\)の\(a=\frac{-1}{CR}\)のパター...
*** 微分方程式(2)の解 [#xa96800d]
次に微分方程式(2)に対しても同様にやってみます。
sageへの入力:
#pre{{
# 微分方程式(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)
}}
&color(blue){laplace(V(t), t, s) == (C*R*s*V(0) + 1)/(C*R...
#pre{{
# 得られた解に初期値V(0)=0を代入し
s22 = s2[0].rhs().subs_expr(V(0) == 0); show(s22)
# 逆ラプラス変換する
inverse_laplace(s22, s, t)
}}
&color(blue){\(\frac{1}{C R s^{2} + s}\)};
&color(blue){-e^(-t/(C*R)) + 1};
微分方程式(1)と同様微分方程式(2)も期待した解を得ることが...
** maximaを使った微分方程式の解法 [#cca4a752]
sageでもある程度の微分方程式は解くことができますが、初期...
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の入力：
#pre{{
# 変数x, yを定義
var('x y')
# maximaを使って微分方程式を定義：diffの前の'に注意
deq = maxima("x^2*'diff(y,x) + 3*x*y = sin(x)/x")
}}
#pre{{
# 微分方程式を解く
d2 = deq.ode2(y,x).ic1(x=pi,y=0)
# 解を表示
show(deq.ode2(y,x))
# maxima の結果をsageの表現に変換するには_sage()_を使い、...
d3 = d2._sage_(); d3.rhs()
}}
&color(blue){\(y={{{\it c}-\cos x}\over{x^3}}\)};
&color(blue){-(cos(x) + 1)/x^3};
#pre{{
# ic1を使って初期値を指定して微分方程式を解く
d4 = deq.ode2(y,x).ic1(x=pi,y=0); d4
}}
&color(blue){y=-(cos(x)+1)/x^3};
#pre{{
# 結果をプロットする
plot(d4._sage_().rhs(), [x, 0, pi])
}}
&ref(sage0-2.png);
このようにsageを使って微分方程式を様々な手法で解くことに...
是非、この機会にsageを使って微分方程式を解いてみてくださ...
** コメント [#k0468c24]
#vote(おもしろかった[10],そうでもない[2],わかりずらい[4])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2009/11/01からのアクセス回数 &counter;
sageを使ってエンジニアにとって身近な微分方程式を解いてみ...
また、実際にsageを使って出力されたワークシートは、以下のU...
http://www.sagenb.org/home/pub/906/
** RC回路 [#f65f28fa]
*** RC直列回路の放電 [#u93cc6d1]
以下のような抵抗（R）とコンデンサー（C）の回路((図と式の...
&ref(RC1.gif);
- 抵抗の両端の電圧をVr
- コンデンサーの両端の電圧をVc
とし、t=0でVcにv0の電圧を掛け、すぐに0にした場合のVc電圧...
キルヒホッフの法則から回路を一周したときの電圧降下は０と...
$$
\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}\...
$$
\begin{array}{l l l}
V(t) + RC\frac{dV(t)}{dt} = 0 & \cdots & (1)
\end{array}
$$
となります。
これをsageを使って解くと以下のようになります。微分方程式...
desolve(微分方程式, [求める関数と変数のリスト], [初期値...
sageへの入力：
#pre{{
# 使用する変数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)
}}
&color(blue){\(e^{-\frac{t}{C R}}\)};
#pre{{
# 解をR=1, C=1のVの変化をプロットする
f(t,R,C) = des
plot(f(t,1,1),[t,0,5])
}}
&ref(sage0.png);
*** RC直列回路の充電 [#z90b8582]
つぎに、t=0から電圧\(V_0\)を掛けて充電した場合のVcの変化...
&ref(RC2.gif);
$$
\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)では右辺が０であり、このような方程式は斉次方程式と呼...
非斉次方程式の一般解 ＝ 斉次方程式の一般解 ＋ 非斉次...
から求めることができます。
- 非斉次方程式の一つの解（特解）は、スイッチを入れて長い...
- desolveの結果から、斉次方程式の一般解は、\(V(t) = A e^{...
t=0でV(t)=0であるとすると、A=−V0 となることが分かります。...
$$
\begin{array}{l l l}
V(t) = V_0(1-e^{\frac{-t}{RC}}) & \cdots & (3)
\end{array}
$$
となります。
sageへの入力:
#pre{{
v(t, R, C) = (1 - des)
plot(v(t, 1, 1),[t,0,5])
}}
&ref(sage0-1.png);
** ラブラス変換を使った微分方程式の解法 [#aa81fa10]
ラプラス変換を使用して微分方程式の解を求めることができま...
*** ラブラス変換とは [#ta611804]
ラプラス変換とは、関数\(f(t)\)に\(e^{-st}\) を掛け、t に...
$$
F(s) = \mathcal{L}(f(t)) = \int_{0}^{\infty}e^{-st}f(t)dt
$$
ここで、\(s\)は\(t\)に無関係な変数であり、\(t\)は実変数で...
sageへの入力:
#pre{{
var('t s')
f = function('f', t)
# ラプラス変換のsageでの表現
view(laplace(f, t, s))
}}
&color(blue){\(\mathcal{L}\left(f\left(t\right), t, s\rig...
*** 一般的なラプラス変換の例 [#u005628c]
主な関数のラプラス変換をした結果を以下に示します。
|元の関数|ラプラス変換結果|h
|\(\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(\omega t)\)|\(\frac{s}{{(\omega^{2}+s^{2})}}\)|
|\(sin(\omega t)\)|\(\frac{\omega}{{(\omega^{2} + s^{2})}...
|\(e^{a t}\)|\(-\frac{1}{{(a - s)}}\)|
|\(t e^{a t}\)|\(\frac{1}{{(a - s)}^{2}}\)|
上記の表を出力するsageの入力：
#pre{{
# 変数の宣言と制約条件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 t^n, laplace(t^n, t, s)
print cos(omega*t), laplace(cos(omega*t), t, s)
print sin(omega*t), laplace(sin(omega*t), t, s)
print exp(a*t), laplace(exp(a*t), t, s)
print t*exp(a*t), laplace(t*exp(a*t), t, s)
}}
出力：
#pre{{
\delta 1
1 1/s
t s^(-2)
t^2 2/s^3
t^n s^(-n - 1)*gamma(n + 1)
cos(omega*t) s/(omega^2 + s^2)
sin(omega*t) omega/(omega^2 + s^2)
e^(a*t) -1/(a - s)
t*e^(a*t) (a - s)^(-2)
}}
*** ラプラス変換の性質 [#bd46fa5d]
ラプラス変換の性質で、特質すべきはラプラス変換の微分によ...
$$
\mathcal{L}(f') = s\mathcal{L}(f) + f(0)
$$
これをsageを使って表現すると以下のようになります。
sage入力：
#pre{{
f = function('f', t)
laplace(diff(f,t), t, s)
}}
&color(blue){s*laplace(f(t), t, s) - f(0)};
*** 微分方程式の解法 [#pdc59ce1]
ラブラス変換を使った微分方程式の解法は、以下の手順で行い...
+ 微分方程式にラプラス変換を施し、微分方程式を変数t領域か...
+ 得られた方程式は演算子$s$の方程式となり、solve関数によ...
+ 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換...
それでは、sageで微分方程式(1)をラプラス変換を使って解くと...
sageの入力：
#pre{{
# deをラプラス変換する
l1 = laplace(de, t, s); l1
}}
&color(blue){(s*laplace(V(t), t, s) - V(0))*C*R + laplace...
#pre{{
# この方程式を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
}}
&color(blue){(laplace(V(t), t, s) == C*R*V(0)/(C*R*s + 1)};
&color(blue){C*R/(C*R*s + 1)};
#pre{{
# 得られた解C*R/(C*R*s + 1)を逆ラプラス変換する
is1 = inverse_laplace(s11, s, t); is1
}}
&color(blue){e^(-t/(C*R))};
これで、ラプラス変換からも\(e^{\frac{-t}{RC}}\)の解を得る...
手計算では、ラプラス変換表を使って逆ラプラス変換を行って...
- \(\frac{CR}{(CRs + 1}\)を\(-CR\)で割ると\(\frac{-1}{(\f...
- ラプラス変換表の\(e^{at}\)の\(a=\frac{-1}{CR}\)のパター...
*** 微分方程式(2)の解 [#xa96800d]
次に微分方程式(2)に対しても同様にやってみます。
sageへの入力:
#pre{{
# 微分方程式(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)
}}
&color(blue){laplace(V(t), t, s) == (C*R*s*V(0) + 1)/(C*R...
#pre{{
# 得られた解に初期値V(0)=0を代入し
s22 = s2[0].rhs().subs_expr(V(0) == 0); show(s22)
# 逆ラプラス変換する
inverse_laplace(s22, s, t)
}}
&color(blue){\(\frac{1}{C R s^{2} + s}\)};
&color(blue){-e^(-t/(C*R)) + 1};
微分方程式(1)と同様微分方程式(2)も期待した解を得ることが...
** maximaを使った微分方程式の解法 [#cca4a752]
sageでもある程度の微分方程式は解くことができますが、初期...
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の入力：
#pre{{
# 変数x, yを定義
var('x y')
# maximaを使って微分方程式を定義：diffの前の'に注意
deq = maxima("x^2*'diff(y,x) + 3*x*y = sin(x)/x")
}}
#pre{{
# 微分方程式を解く
d2 = deq.ode2(y,x).ic1(x=pi,y=0)
# 解を表示
show(deq.ode2(y,x))
# maxima の結果をsageの表現に変換するには_sage()_を使い、...
d3 = d2._sage_(); d3.rhs()
}}
&color(blue){\(y={{{\it c}-\cos x}\over{x^3}}\)};
&color(blue){-(cos(x) + 1)/x^3};
#pre{{
# ic1を使って初期値を指定して微分方程式を解く
d4 = deq.ode2(y,x).ic1(x=pi,y=0); d4
}}
&color(blue){y=-(cos(x)+1)/x^3};
#pre{{
# 結果をプロットする
plot(d4._sage_().rhs(), [x, 0, pi])
}}
&ref(sage0-2.png);
このようにsageを使って微分方程式を様々な手法で解くことに...
是非、この機会にsageを使って微分方程式を解いてみてくださ...
** コメント [#k0468c24]
#vote(おもしろかった[10],そうでもない[2],わかりずらい[4])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
