sage/graph/Interface201511(サーミスタ温度特性のグラフ)
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2015/10/03からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/66/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう:Interface 2015/11 [#wf062...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
Interface 2015/11号にMathematicaを使った例題が紹介されて...
これをSageを使って試してみます。
最初に図7(Interface 2015/11から引用)に紹介されている例...
&ref(fig-7.png);
** 10ビット三角関数テーブルを作る [#w04e542c]
図7の最初の部分は、sinのテーブルを作成することです。
Sageのリストでは、リストの要素単位に演算をできません。
このような場合には、numpyのarrayを使います。
pythonでリストを作るときには、以下の書式のリスト内包を使...
#pre{{
[ 式 for 変数 in リスト ]
}}
forで変数に値をセットするためのリストには、通常rangeやara...
今回は最後の2πを含めたリストを作るために、0から2πを均等分...
sageへの入力:
#pre{{
# 要素単位への演算をするため、numpyのarrayを使用
# arangeだと最後の0が含まれないのでlinspaceで21個出力する...
import numpy as np
np.array([sin(N(x)) for x in np.linspace(0, 2*pi, 21)])
}}
#pre{{
array([ 0.00000000e+00, 3.09016994e-01, 5.87785252e-...
8.09016994e-01, 9.51056516e-01, 1.00000000e+...
9.51056516e-01, 8.09016994e-01, 5.87785252e-...
3.09016994e-01, 1.22464680e-16, -3.09016994e-...
-5.87785252e-01, -8.09016994e-01, -9.51056516e-...
-1.00000000e+00, -9.51056516e-01, -8.09016994e-...
-5.87785252e-01, -3.09016994e-01, -2.44929360e-...
}}
直前の結果を参照するには、_(アンダーバー)を使います。
以下の様に直前の結果を1000倍し、round関数で小数点以下を取...
10ビット三角関数テーブルを作ります。
sageへの入力:
#pre{{
# 直前の結果は、_で参照できる
1023*_
}}
#pre{{
array([ 0.00000000e+00, 3.16124385e+02, 6.01304313e+...
8.27624385e+02, 9.72930816e+02, 1.02300000e+...
9.72930816e+02, 8.27624385e+02, 6.01304313e+...
3.16124385e+02, 1.25281368e-13, -3.16124385e+...
-6.01304313e+02, -8.27624385e+02, -9.72930816e+...
-1.02300000e+03, -9.72930816e+02, -8.27624385e+...
-6.01304313e+02, -3.16124385e+02, -2.50562735e-...
}}
sageへの入力:
#pre{{
round(_)
}}
#pre{{
array([ 0., 316., 601., 828., 973., 1023., ...
601., 316., 0., -316., -601., -828., -...
-973., -828., -601., -316., -0.])
}}
*** 結果をグラフにプロットする [#yf727f40]
リストをグラフにプロットするには、list_plot関数を使います。
そのままだとグラフが大きく表示されるので、figsize=5で調節...
sageへの入力:
#pre{{
list_plot(_, figsize=5)
}}
&ref(sage0.png);
最初のリスト内包の式で、sin(N(x))としていましたが、sin(x)...
Sageでは、数式がそのまま使われるので、分数やsin等の関数が...
sageへの入力:
#pre{{
# Nを使わないと、数式のまま表示される
[sin(x) for x in np.arange(0, 2*pi, pi/10)]
}}
#pre{{
[0,
1/4*sqrt(5) - 1/4,
sin(1/5*pi),
1/4*sqrt(5) + 1/4,
sin(2/5*pi),
1,
sin(3/5*pi),
1/4*sqrt(5) + 1/4,
sin(4/5*pi),
1/4*sqrt(5) - 1/4,
0,
-1/4*sqrt(5) + 1/4,
sin(6/5*pi),
-1/4*sqrt(5) - 1/4,
sin(7/5*pi),
-1,
sin(8/5*pi),
-1/4*sqrt(5) - 1/4,
sin(9/5*pi),
-1/4*sqrt(5) + 1/4]
}}
リストの部分列を取り出すには、[始まりのインデクス:終わり...
以下は、aから最初の0, 1, 2番目の要素を取り出すために、a[0...
丸括弧()で括られている部分は、タプルと呼ばれ、座標点など...
sageへの入力:
#pre{{
# pythonは0始まりで、終わりのインデックスの前まで出力する...
a = [(N(x), sin(N(x))) for x in np.linspace(0, 2*pi, 21)]...
}}
#pre{{
[(0.000000000000000, 0.000000000000000),
(0.314159265358979, 0.309016994374947),
(0.628318530717959, 0.587785252292473)]
}}
*** データフィッティング [#qdcc8058]
Sageでデータに最も合ったモデルを求める関数がfind_fit関数...
find_fitには、データとモデルを入力としますが、モデルに使...
以下では5次関数をモデルとする、model(x)関数を変数(係数)...
solution_dict=Trueは、結果を辞書型で返すように指定します...
find_fitの結果を使った関数がとても簡単に求まります。
sageへの入力:
#pre{{
(x, a0, a1, a2, a3, a4, a5) = var('x a0 a1 a2 a3 a4 a5')
model(x) = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5
fit = find_fit(a, model, solution_dict=True)
#show(model)
#show(fit)
print model
print fit
}}
*** グラフの重ね合わせ [#p317bf59]
データのリストプロットとフィッティングで求めた関数b(x)を...
プロット関数の結果を+記号で足し合わせます。
関数のプロットには、plot関数を使い、[x, 0, 2*pi]でプロッ...
rgbcolor='green'で線色を緑に指定します。
図7ではSin関数とフィッティング曲線の重ね合わせですが、こ...
sageへの入力:
#pre{{
b(x) = model.subs(fit)
list_plot(a, figsize=5)+plot(b(x), [x, 0, 2*pi], rgbcolor...
}}
&ref(sage1.png);
次に、Sin曲線とフィッティング曲線の差をプロットします。pl...
b(x)-sin(x)をそのまま書くとSin曲線とフィッティング曲線の...
sageへの入力:
#pre{{
plot(b(x)-sin(x), [x, 0, 2*pi], figsize=5)
}}
&ref(sage2.png);
最初の三角関数テーブルを1行のリスト内包で書くと以下の様に...
sageへの入力:
#pre{{
# pythonのリスト内包の表現は強力です
[N(round(1023*sin(x))) for x in np.linspace(0, 2*pi, 21)]
}}
#pre{{
[0.000000000000000,
316.000000000000,
601.000000000000,
828.000000000000,
973.000000000000,
1023.00000000000,
973.000000000000,
828.000000000000,
601.000000000000,
316.000000000000,
0.000000000000000,
-316.000000000000,
-601.000000000000,
-828.000000000000,
-973.000000000000,
-1023.00000000000,
-973.000000000000,
-828.000000000000,
-601.000000000000,
-316.000000000000,
0.000000000000000]
}}
** 音声データ [#xa17ed32]
図8(Interface 2015/11から引用)では、Mathematicaを使った...
&ref(fig-8.png);
Sageには音声データを扱う関数がないので、Pythonのライブラ...
今回は、音声データの取り込みにSoundFileライブラリを使用し...
波形とスペクトログラムの表示を紹介するのに留めます。
*** 準備 [#a79036b5]
最初に、必要なライブラリをインポートします。
sageへの入力:
#pre{{
# 音声データの処理用パッケージ
import soundfile as sf
import pylab as pl
import scipy.fftpack as Fourier
from scipy.signal import butter, lfilter, freqz
import time
}}
波形とスペクトログラムをプロットする関数と音声データの書...
sageへの入力:
#pre{{
# pylabを使ったスペックグラムと波形の処理用に関数を定義
# 画像ファイルの保存と表示関数
def saveAndShowPlot(filename, fac=0.75):
output = DATA + filename
pl.savefig(output)
width = int(640*fac)
html('<img src="%s?%s" width="%spx">'%(filename, time...
# specgramと波形のプロット関数
def plotSpecgram(wave, Fs):
# FFTのサンプル数
N = 512
# FFTで用いるハミング窓
hammingWindow = np.hamming(N)
# specgramのプロット
t = range(0, len(wave))
pl.figure()
pl.subplot(2, 1, 1)
pl.specgram(wave, NFFT=N, Fs=Fs, noverlap=N/2, window...
# waveのプロット
pl.subplot(2, 1, 2)
pl.plot(t, wave, 'b')
pl.grid()
# ファイル保存
def recordSound(sources, name, Hz):
filename = DATA + name
sf.write(sources, filename, samplerate=Hz)
# 結果を確認する方法
html('<a href="%s">録音結果(右クリックでファイルをダ...
}}
*** サウンドデータの作成 [#u14b218d]
サウンドデータとして1秒間で20Hz〜20kHzに上昇するスイープ...
これもPythonのリスト内包を使って簡単にできます。
sageへの入力:
#pre{{
smpRate = 44100 # サンプリングレート44.1kHz
fLo = 20
fHi = 20000
wave = [N(sin(2*pi*(fLo+(fHi/2-fLo)*t)*t)) for t in np.li...
}}
音声データのプロットと保存用に定義したplotSpecgram関数とs...
生成したサウンドデータの波形とスペクトログラムを表示し、...
sageへの入力:
#pre{{
plotSpecgram(wave, smpRate)
saveAndShowPlot('wave.png')
}}
&ref(th_wave.jpg);
sageへの入力:
#pre{{
# waveをファイルに保存
recordSound(wave, 'swep441k20-20k.wav', smpRate)
}}
- 録音結果(&ref(swep441k20-20k.wav);)
スペクトログラムの形が図8(ソノグラフ)と異なるので、以...
のスペクトログラムをプロットしてみました。正しく表示され...
の定義が異なるのかも知れませんが、ここではこれ以上深入り...
sageへの入力:
#pre{{
# specgramの形が異なるため、以前Interfaceで紹介されていた...
# ファイルからの読み込み
rec, Fs = sf.read(DATA+"voice_howling.wav")
sample =rec.flatten()
}}
sageへの入力:
#pre{{
# 以前と同じ結果がでたので、問題なしとする
plotSpecgram(sample, Fs)
saveAndShowPlot('sample.png')
}}
&ref(th_sample.jpg);
** 連立方程式を使った応用例 [#eeb82567]
図11(Interface 2015/11から引用)では、サーミスタの抵抗と...
から\(v_0\)と\(T_s\)の関係式を求めています。
$$
\begin{eqnarray}
eq_1 & = & R_s == R_2 e^{B \left ( \frac{1}{T_s} - \frac...
eq_2 & = & v_0 == \frac{R}{R_s + R} v_e
\end{eqnarray}
$$
最初に使用する変数R, Rs, R2, B, Ts, T2, ve, v0を変数定義...
sageへの入力:
#pre{{
R, Rs, R2, B, Ts, T2, ve, v0 = var('R Rs R2 B Ts T2 ve v0')
# 以下の関係式が与えられた場合のTsとv0の関係を求める
eq1 = Rs == R2*exp(B*(1/Ts - 1/T2))
eq2 = v0 == R/(Rs + R)*ve
}}
*** solve関数の使い方 [#w377bc39]
solve関数を使って連立方程式を解きます。
本来なら、以下の様にすれば解けるのですが、
#pre{{
solve([eq1, eq2], [Ts,Rs])
}}
eq2にTsが含まれていないため、Sageでは上手く解くことができ...
ここでは、連立方程式の解法を順番にSageを使って解を求めて...
最初にeq2からRsを求め、その結果をRs_eqにセットします。
sageへの入力:
#pre{{
# eq2をRsについて整理
Rs_eq = solve(eq2, Rs); Rs_eq
}}
#pre{{
[Rs == -R*(v0 - ve)/v0]
}}
つぎにRs_eqをeq1に代入し、Tsを求め、右辺の式を取りだしま...
solve関数の結果はリスト形式で返るので、[0]で最初の1個を取...
求まった解は、図11の表現とは異なりますが、logに-1を掛ける...
$$
\begin{eqnarray}
T_s & = & -\frac{B T_2}{T_2 log \left [ -\frac{R_2 v_0}{R...
& = & -\frac{B T_2}{T_2 log \left [ \frac{R_2 v_0}{R v_e ...
& = & \frac{-B T_2}{-T_2 log \left [ \frac{R v_e - R v_0 ...
& = & \frac{B T_2}{B + T_2 log \left [ \frac{R (v_e - v_0...
\end{eqnarray}
$$
sageへの入力:
#pre{{
# Rs_eqの結果で、eq1を解く
sol = solve(eq1.subs(Rs_eq[0]), Ts)[0].rhs()
#show(sol)
sol
}}
#pre{{
-B*T2/(T2*log(-R2*v0/(R*v0 - R*ve)) - B)
}}
次に変数R, ve, R2, T2, Bに値をセットし、その時のTsとv0の...
これには、変数名:値の対をカンマで区切り、{}で囲んだ辞書型...
これで目的とするTsを求める関数tsFunc(v0)が求まりました。
plot関数を使ってtsFuncをプロットすると図11のTsとv0の関係...
sageへの入力:
#pre{{
# 他の変数に値をセットして、Tsを変数v0で定義した関数にする
tsFunc(v0) = sol.subs({R:10000, ve:1, R2: 30000, T2: 273,...
#show(tsFunc)
tsFunc
}}
#pre{{
v0 |--> -311220/(91*log(-3*v0/(v0 - 1)) - 1140)
}}
sageへの入力:
#pre{{
plot(tsFunc(v0)-273, [v0, 0.1, 0.9], figsize=5)
}}
&ref(sage3.png);
** コメント [#mb2dcef5]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2015/10/03からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/66/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう:Interface 2015/11 [#wf062...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
Interface 2015/11号にMathematicaを使った例題が紹介されて...
これをSageを使って試してみます。
最初に図7(Interface 2015/11から引用)に紹介されている例...
&ref(fig-7.png);
** 10ビット三角関数テーブルを作る [#w04e542c]
図7の最初の部分は、sinのテーブルを作成することです。
Sageのリストでは、リストの要素単位に演算をできません。
このような場合には、numpyのarrayを使います。
pythonでリストを作るときには、以下の書式のリスト内包を使...
#pre{{
[ 式 for 変数 in リスト ]
}}
forで変数に値をセットするためのリストには、通常rangeやara...
今回は最後の2πを含めたリストを作るために、0から2πを均等分...
sageへの入力:
#pre{{
# 要素単位への演算をするため、numpyのarrayを使用
# arangeだと最後の0が含まれないのでlinspaceで21個出力する...
import numpy as np
np.array([sin(N(x)) for x in np.linspace(0, 2*pi, 21)])
}}
#pre{{
array([ 0.00000000e+00, 3.09016994e-01, 5.87785252e-...
8.09016994e-01, 9.51056516e-01, 1.00000000e+...
9.51056516e-01, 8.09016994e-01, 5.87785252e-...
3.09016994e-01, 1.22464680e-16, -3.09016994e-...
-5.87785252e-01, -8.09016994e-01, -9.51056516e-...
-1.00000000e+00, -9.51056516e-01, -8.09016994e-...
-5.87785252e-01, -3.09016994e-01, -2.44929360e-...
}}
直前の結果を参照するには、_(アンダーバー)を使います。
以下の様に直前の結果を1000倍し、round関数で小数点以下を取...
10ビット三角関数テーブルを作ります。
sageへの入力:
#pre{{
# 直前の結果は、_で参照できる
1023*_
}}
#pre{{
array([ 0.00000000e+00, 3.16124385e+02, 6.01304313e+...
8.27624385e+02, 9.72930816e+02, 1.02300000e+...
9.72930816e+02, 8.27624385e+02, 6.01304313e+...
3.16124385e+02, 1.25281368e-13, -3.16124385e+...
-6.01304313e+02, -8.27624385e+02, -9.72930816e+...
-1.02300000e+03, -9.72930816e+02, -8.27624385e+...
-6.01304313e+02, -3.16124385e+02, -2.50562735e-...
}}
sageへの入力:
#pre{{
round(_)
}}
#pre{{
array([ 0., 316., 601., 828., 973., 1023., ...
601., 316., 0., -316., -601., -828., -...
-973., -828., -601., -316., -0.])
}}
*** 結果をグラフにプロットする [#yf727f40]
リストをグラフにプロットするには、list_plot関数を使います。
そのままだとグラフが大きく表示されるので、figsize=5で調節...
sageへの入力:
#pre{{
list_plot(_, figsize=5)
}}
&ref(sage0.png);
最初のリスト内包の式で、sin(N(x))としていましたが、sin(x)...
Sageでは、数式がそのまま使われるので、分数やsin等の関数が...
sageへの入力:
#pre{{
# Nを使わないと、数式のまま表示される
[sin(x) for x in np.arange(0, 2*pi, pi/10)]
}}
#pre{{
[0,
1/4*sqrt(5) - 1/4,
sin(1/5*pi),
1/4*sqrt(5) + 1/4,
sin(2/5*pi),
1,
sin(3/5*pi),
1/4*sqrt(5) + 1/4,
sin(4/5*pi),
1/4*sqrt(5) - 1/4,
0,
-1/4*sqrt(5) + 1/4,
sin(6/5*pi),
-1/4*sqrt(5) - 1/4,
sin(7/5*pi),
-1,
sin(8/5*pi),
-1/4*sqrt(5) - 1/4,
sin(9/5*pi),
-1/4*sqrt(5) + 1/4]
}}
リストの部分列を取り出すには、[始まりのインデクス:終わり...
以下は、aから最初の0, 1, 2番目の要素を取り出すために、a[0...
丸括弧()で括られている部分は、タプルと呼ばれ、座標点など...
sageへの入力:
#pre{{
# pythonは0始まりで、終わりのインデックスの前まで出力する...
a = [(N(x), sin(N(x))) for x in np.linspace(0, 2*pi, 21)]...
}}
#pre{{
[(0.000000000000000, 0.000000000000000),
(0.314159265358979, 0.309016994374947),
(0.628318530717959, 0.587785252292473)]
}}
*** データフィッティング [#qdcc8058]
Sageでデータに最も合ったモデルを求める関数がfind_fit関数...
find_fitには、データとモデルを入力としますが、モデルに使...
以下では5次関数をモデルとする、model(x)関数を変数(係数)...
solution_dict=Trueは、結果を辞書型で返すように指定します...
find_fitの結果を使った関数がとても簡単に求まります。
sageへの入力:
#pre{{
(x, a0, a1, a2, a3, a4, a5) = var('x a0 a1 a2 a3 a4 a5')
model(x) = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5
fit = find_fit(a, model, solution_dict=True)
#show(model)
#show(fit)
print model
print fit
}}
*** グラフの重ね合わせ [#p317bf59]
データのリストプロットとフィッティングで求めた関数b(x)を...
プロット関数の結果を+記号で足し合わせます。
関数のプロットには、plot関数を使い、[x, 0, 2*pi]でプロッ...
rgbcolor='green'で線色を緑に指定します。
図7ではSin関数とフィッティング曲線の重ね合わせですが、こ...
sageへの入力:
#pre{{
b(x) = model.subs(fit)
list_plot(a, figsize=5)+plot(b(x), [x, 0, 2*pi], rgbcolor...
}}
&ref(sage1.png);
次に、Sin曲線とフィッティング曲線の差をプロットします。pl...
b(x)-sin(x)をそのまま書くとSin曲線とフィッティング曲線の...
sageへの入力:
#pre{{
plot(b(x)-sin(x), [x, 0, 2*pi], figsize=5)
}}
&ref(sage2.png);
最初の三角関数テーブルを1行のリスト内包で書くと以下の様に...
sageへの入力:
#pre{{
# pythonのリスト内包の表現は強力です
[N(round(1023*sin(x))) for x in np.linspace(0, 2*pi, 21)]
}}
#pre{{
[0.000000000000000,
316.000000000000,
601.000000000000,
828.000000000000,
973.000000000000,
1023.00000000000,
973.000000000000,
828.000000000000,
601.000000000000,
316.000000000000,
0.000000000000000,
-316.000000000000,
-601.000000000000,
-828.000000000000,
-973.000000000000,
-1023.00000000000,
-973.000000000000,
-828.000000000000,
-601.000000000000,
-316.000000000000,
0.000000000000000]
}}
** 音声データ [#xa17ed32]
図8(Interface 2015/11から引用)では、Mathematicaを使った...
&ref(fig-8.png);
Sageには音声データを扱う関数がないので、Pythonのライブラ...
今回は、音声データの取り込みにSoundFileライブラリを使用し...
波形とスペクトログラムの表示を紹介するのに留めます。
*** 準備 [#a79036b5]
最初に、必要なライブラリをインポートします。
sageへの入力:
#pre{{
# 音声データの処理用パッケージ
import soundfile as sf
import pylab as pl
import scipy.fftpack as Fourier
from scipy.signal import butter, lfilter, freqz
import time
}}
波形とスペクトログラムをプロットする関数と音声データの書...
sageへの入力:
#pre{{
# pylabを使ったスペックグラムと波形の処理用に関数を定義
# 画像ファイルの保存と表示関数
def saveAndShowPlot(filename, fac=0.75):
output = DATA + filename
pl.savefig(output)
width = int(640*fac)
html('<img src="%s?%s" width="%spx">'%(filename, time...
# specgramと波形のプロット関数
def plotSpecgram(wave, Fs):
# FFTのサンプル数
N = 512
# FFTで用いるハミング窓
hammingWindow = np.hamming(N)
# specgramのプロット
t = range(0, len(wave))
pl.figure()
pl.subplot(2, 1, 1)
pl.specgram(wave, NFFT=N, Fs=Fs, noverlap=N/2, window...
# waveのプロット
pl.subplot(2, 1, 2)
pl.plot(t, wave, 'b')
pl.grid()
# ファイル保存
def recordSound(sources, name, Hz):
filename = DATA + name
sf.write(sources, filename, samplerate=Hz)
# 結果を確認する方法
html('<a href="%s">録音結果(右クリックでファイルをダ...
}}
*** サウンドデータの作成 [#u14b218d]
サウンドデータとして1秒間で20Hz〜20kHzに上昇するスイープ...
これもPythonのリスト内包を使って簡単にできます。
sageへの入力:
#pre{{
smpRate = 44100 # サンプリングレート44.1kHz
fLo = 20
fHi = 20000
wave = [N(sin(2*pi*(fLo+(fHi/2-fLo)*t)*t)) for t in np.li...
}}
音声データのプロットと保存用に定義したplotSpecgram関数とs...
生成したサウンドデータの波形とスペクトログラムを表示し、...
sageへの入力:
#pre{{
plotSpecgram(wave, smpRate)
saveAndShowPlot('wave.png')
}}
&ref(th_wave.jpg);
sageへの入力:
#pre{{
# waveをファイルに保存
recordSound(wave, 'swep441k20-20k.wav', smpRate)
}}
- 録音結果(&ref(swep441k20-20k.wav);)
スペクトログラムの形が図8(ソノグラフ)と異なるので、以...
のスペクトログラムをプロットしてみました。正しく表示され...
の定義が異なるのかも知れませんが、ここではこれ以上深入り...
sageへの入力:
#pre{{
# specgramの形が異なるため、以前Interfaceで紹介されていた...
# ファイルからの読み込み
rec, Fs = sf.read(DATA+"voice_howling.wav")
sample =rec.flatten()
}}
sageへの入力:
#pre{{
# 以前と同じ結果がでたので、問題なしとする
plotSpecgram(sample, Fs)
saveAndShowPlot('sample.png')
}}
&ref(th_sample.jpg);
** 連立方程式を使った応用例 [#eeb82567]
図11(Interface 2015/11から引用)では、サーミスタの抵抗と...
から\(v_0\)と\(T_s\)の関係式を求めています。
$$
\begin{eqnarray}
eq_1 & = & R_s == R_2 e^{B \left ( \frac{1}{T_s} - \frac...
eq_2 & = & v_0 == \frac{R}{R_s + R} v_e
\end{eqnarray}
$$
最初に使用する変数R, Rs, R2, B, Ts, T2, ve, v0を変数定義...
sageへの入力:
#pre{{
R, Rs, R2, B, Ts, T2, ve, v0 = var('R Rs R2 B Ts T2 ve v0')
# 以下の関係式が与えられた場合のTsとv0の関係を求める
eq1 = Rs == R2*exp(B*(1/Ts - 1/T2))
eq2 = v0 == R/(Rs + R)*ve
}}
*** solve関数の使い方 [#w377bc39]
solve関数を使って連立方程式を解きます。
本来なら、以下の様にすれば解けるのですが、
#pre{{
solve([eq1, eq2], [Ts,Rs])
}}
eq2にTsが含まれていないため、Sageでは上手く解くことができ...
ここでは、連立方程式の解法を順番にSageを使って解を求めて...
最初にeq2からRsを求め、その結果をRs_eqにセットします。
sageへの入力:
#pre{{
# eq2をRsについて整理
Rs_eq = solve(eq2, Rs); Rs_eq
}}
#pre{{
[Rs == -R*(v0 - ve)/v0]
}}
つぎにRs_eqをeq1に代入し、Tsを求め、右辺の式を取りだしま...
solve関数の結果はリスト形式で返るので、[0]で最初の1個を取...
求まった解は、図11の表現とは異なりますが、logに-1を掛ける...
$$
\begin{eqnarray}
T_s & = & -\frac{B T_2}{T_2 log \left [ -\frac{R_2 v_0}{R...
& = & -\frac{B T_2}{T_2 log \left [ \frac{R_2 v_0}{R v_e ...
& = & \frac{-B T_2}{-T_2 log \left [ \frac{R v_e - R v_0 ...
& = & \frac{B T_2}{B + T_2 log \left [ \frac{R (v_e - v_0...
\end{eqnarray}
$$
sageへの入力:
#pre{{
# Rs_eqの結果で、eq1を解く
sol = solve(eq1.subs(Rs_eq[0]), Ts)[0].rhs()
#show(sol)
sol
}}
#pre{{
-B*T2/(T2*log(-R2*v0/(R*v0 - R*ve)) - B)
}}
次に変数R, ve, R2, T2, Bに値をセットし、その時のTsとv0の...
これには、変数名:値の対をカンマで区切り、{}で囲んだ辞書型...
これで目的とするTsを求める関数tsFunc(v0)が求まりました。
plot関数を使ってtsFuncをプロットすると図11のTsとv0の関係...
sageへの入力:
#pre{{
# 他の変数に値をセットして、Tsを変数v0で定義した関数にする
tsFunc(v0) = sol.subs({R:10000, ve:1, R2: 30000, T2: 273,...
#show(tsFunc)
tsFunc
}}
#pre{{
v0 |--> -311220/(91*log(-3*v0/(v0 - 1)) - 1140)
}}
sageへの入力:
#pre{{
plot(tsFunc(v0)-273, [v0, 0.1, 0.9], figsize=5)
}}
&ref(sage3.png);
** コメント [#mb2dcef5]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc