スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Scilabでおもりの吊り下げ

順番に行くならば、微分方程式による物理現象のモデル化(PDF)の次の問題は3個の質点の吊り下げ問題です。

001_20130729204658.png
図35 3個のおもりの吊り下げ


図35に示したように,3個のおもりが,4本のいとでつながれ,天井からつり下げられている.この時の,4本の糸にかかる張力と糸が水平面となす角度Ti(i=1,2,3,4)を求めなさい.


これは、いままでのような常微分方程式を解くものではなく、非線形連立方程式の数値解を求める問題です。
微分方程式による物理現象のモデル化(PDF)では自分で非線形方程式ソルバのプログラムをしていますが、ScilabやOctaveには非線形連立方程式の数値解を求めるための命令fsolveが実装されているので、実は簡単に数値解を求めることが出来ます。

著者の方もこのことに気づいたようでOctaveの精義では、ほとんど同じ問題がfsolveで計算されています。
したがって、今回はOctaveの精義のほうに倣ってfsolveで解くことにします。

今回は、この問題をScilabのfsolveを利用して解く事を考えます。
さらに、おもりの数が増えたり減ったりしたときも計算できるようにn個のおもりをつるしたときの計算が出来るようにします。


非線形連立方程式


解くべき方程式は力のつりあい、x方向とy方向の距離から以下連立方程式が得られます。

T_i \cos\theta_i - T_{i+1} \cos\theta_{i+1} = 0

T_i \sin\theta_i - T_{i+1} \sin\theta_{i+1} - W_i = 0

\sum_{i=1}^n L_i \cos\theta_i - L = 0

\sum_{i=1}^n L_i \sin\theta_i = 0

おもりの数がn個のとき、方程式は2n+2個、求める変数も張力Tiがn+1個と角度θiがn+1個でやはり計2n+2個になります。

プログラムを書く上では「2種類の変数がn+1個ずつ(つまりn+1個の要素を持つベクトル2つ)」というよりも「2n+2個の要素を持つベクトル1つ」のほうが都合が良いのでθiとTiをまとめて実数ベクトルZで以下のようにあらわすことにします。

\theta_i = Z_i

T_i = Z_{(n+1)+i}

行列の形に書き下すと

Z = \begin{pmatrix}Z_1 \\Z_2 \\\vdots \\Z_{n+1} \\Z_{n+2} \\Z_{n+3} \\\vdots \\Z_{2n+2}\end{pmatrix}=\begin{pmatrix}\theta_1 \\\theta_2 \\\vdots \\\theta_{n+1} \\T_{1} \\T_{2} \\\vdots \\T_{n+1}\end{pmatrix}

プログラムを書くときにミスしがちなので、解くべき方程式も丁寧に置き換えて書き直します。

Z_{(n+1)+i}\cos Z_i - Z_{(n+1)+(i+1)}\cos Z_{i+1} = 0

Z_{(n+1)+i}\sin Z_i - Z_{(n+1)+(i+1)}\sin Z_{i+1} - W_i = 0

\sum_{i=1}^{n}L_i\cos Z_i - L = 0

\sum_{i=1}^{n}L_i\sin Z_i = 0

数値計算


以上を踏まえて作成したプログラムがmass_sce.txtです。

002_20130729204658.png

Fig.2: 計算結果。糸の長さはL1=L2=L3=1,L4=3,L=4。おもりの重さはW1=W2=1, W2=2


clear;

// *** 入力パラメータ ***
// おもりの数
n = 3;

// おもりの重さ
for i = 1:n do
W(i) = input(strcat(["W",string(i)," = "]));
end
// 水平方向の距離
l = input("L = ");
// 糸の長さ
for i = 1:(n + 1) do
L(i) = input(strcat(["L",string(i)," = "]));
end

// *** 解くべき連立方程式 ***
function R = f(Z)
R(1:n) = Z(n+2:2*n+1) .* cos(Z(1:n)) - Z(n+3:2*n+2) .* cos(Z(2:n+1));
R(n+1:2*n) = Z(n+2:2*n+1) .* sin(Z(1:n)) - Z(n+3:2*n+2) .* sin(Z(2:n+1)) - W(1:n);
R(2*n+1) = sum(L .* cos(Z(1:(n+1)))) - l;
R(2*n+2) = sum(L .* sin(Z(1:(n+1))));
endfunction
// 初期値
Q0 = ones(2*n+2,1);
Q0(n+1) = -1;
// 非線形方程式ソルバ
Q = fsolve(Q0,f);

// *** おもりの位置の計算とプロット ***
DX = L .* cos(Q(1:n+1));
DY= - L .* sin(Q(1:n+1));
X(1) = DX(1);
Y(1) = DY(1);
for i = 2:n do
X(i) = X(i-1) + DX(i);
Y(i) = Y(i-1) + DY(i);
end
// 縦横比を等しくする
h = scf(); // ウィンドウを作成
ha = h.children(1); // Axes(座標軸)オブジェクトへのハンドルを取得
ha.isoview = "on"; // 座標軸の縦横比を等しくする
ha.data_bounds = [0,-1 * ceil(max([abs(Y);l])); ceil(max([abs(Y);l])),0]; // 座標軸表示範囲の設定
// データのプロット
plot([0;X;l],[0;Y;0]);
plot(X,Y,'or','markersize',10);

180*(Q(1:n+1))/%pi


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 非線形方程式ソルバ 

Scilabで金属の化学ポテンシャル

Scilabの非線形方程式ソルバfsolveを用いて、積分を含む方程式である金属の化学ポテンシャルを計算しました。
対象とする系は、自由電子近似が良く成り立つはずのアルミニウムとし、数値計算による化学ポテンシャルの計算結果とゾンマーフェルト展開を用いた近似値の比較を行いました。

その結果、この系に関しては、数値計算とゾンマーフェルト展開を用いた計算結果がアルミニウムの融点以上までほとんど同じ結果を示すことがわかりました。

001_20130609181459.pngFig.1: アルミニウムの化学ポテンシャルの温度依存性、ゾンマーフェルト展開を用いた近似(青実線)とfsolveを利用した数値計算(赤点線)。インセットは自由電子近似による電子の状態密度(青)、低温時(10K,緑)と高温時(赤)の電子の分布は絶対零度の状態密度に各温度のフェルミ分布関数を掛けたもの、エネルギーの原点はフェルミ準位。



金属電子論と数値計算


Scilabで差し込みグラフ:金属の比熱では実験的に求められたアルミニウムの電子比熱係数γ=1.35(mJ/mol/K^2)という値をもとに電子比熱の温度依存性をグラフにプロットしました。

これに対して、次は理論的にアルミニウムの電子比熱を予想してみます。
差し当たり今回は、電子比熱を計算する前段階として、アルミニウムの化学ポテンシャルを自由電子近似と数値計算から求めます。

金属の電子構造を計算するに当たり、次の2つの近似がよく使われます。
  • 自由電子近似(free electron approximation)
  • ゾンマーフェルト展開(Sommerfeld expansion)

自由電子近似はナトリウムやアルミニウムといった比較的単純な金属に対しては良い近似になることが知られています。また、ゾンマーフェルト展開は温度が充分低い条件で成り立ちます。

ゾンマーフェルト展開は、複雑な複雑な積分の計算を、比較的計算しやすい解析的な式に変形するための級数展開なのですが、Scilabなどで数値計算を行うことでこの近似に頼らずに電子比熱が計算できます。
今回は「自由電子近似+ゾンマーフェルト展開」と「自由電子近似+数値計算」の両方を行い比較します。

金属の化学ポテンシャル


金属の電子比熱は、以下の式で定義されます。

c_e = \frac{\mathrm{d} u_e}{\mathrm{d} T}

ここで、電子系の内部エネルギーは

u_e = \int^{\infty}_{-\infty}\epsilon f(\epsilon) D(\epsilon)\rm{d}\epsilon

さらにD(ε)は電子の状態密度(density of states)でf(ε)はフェルミ・ディラック分布関数です。

f(\epsilon) = \frac{1}{\exp{\left(\frac{\epsilon-\mu}{k_B - T}\right)}+1}

まずD(ε)ですが、これは金属の個性を現す重要なパラメータなのですが、自由電子近似と金属の体積、価電子数から求めることが出来ます(後述)。
一方で、フェルミ・ディラック分布関数内の化学ポテンシャルμは、温度の関数となり少々厄介です。本エントリの主題は、この化学ポテンシャルをScilabの非線形方程式ソルバーを用いて数値的に計算し、ゾンマーフェルト展開によって近似的に求められた化学ポテンシャルと比較を行うことです。

自由電子近似


詳細には立ち入りませんが、自由電子近似による金属の状態密度とフェルミエネルギー(絶対零度における化学ポテンシャル)は、リュードベリ原子単位系では以下の様に表されます。

D(\epsilon) = \frac{V}{2 \pi^2} \epsilon^{\frac{1}{2}}

\epsilon_F = \left( \frac{3 \pi^2 n}{V} \right)^{\frac{2}{3}}

ここでVは単位物質量あたりの体積です。具体的には格子体積やモル体積を用います。今回は差し当たり単位原子辺りの体積(atomic volume)を用いて計算します。nは電子の数で、今回は単位原子辺りの電子の数なので価電子数と一致します。(参考:金持 徹著固体電子論)

自由電子近似は、ナトリウムやアルミニウムでは良く成り立ちますが、遷移金属などには適用できません。

ゾンマーフェルト展開


これも詳細には立ち入らず結論だけ書きます。
ゾンマーフェルト展開の結果、内部エネルギーue、化学ポテンシャルμ、電子比熱ceは

u_e = u_{e0} + \frac{\pi^2}{6}k_B^2 D(\epsilon_F)T^2

\mu = \epsilon_F - \frac{\pi^2}{6}k_B^2 \frac{D^{'}(\epsilon_F)}{D(\epsilon_F)}T^2

c_e = \frac{\pi^2}{3}k_B^2 D(\epsilon_F)T

となります。(参考:無機化学特論I ゾンマーフェルト展開について)

ゾンマーフェルト展開は充分低温のみで適用できます。(参考:Sommerfeld 展開の導出- 武智公平)

自由電子近似+ゾンマーフェルト展開


ゾンマーフェルト展開に自由電子近似の結果を代入すると、金属の化学ポテンシャルや電子比熱の値が具体的に計算できます。

\mu = \epsilon_F \frac{\pi^2 k_B}{12 \epsilon_F}T^2

c_e = \frac{\pi^2 k_B^2 n}{2 \epsilon_F}T

押さえておかなくてはならない重要な点は、「自由電子近似」と「ゾンマーフェルト展開」は相互に全く関係の無い別の近似であり、適用できる条件もそれぞれ異なるということです。

数値計算


さて、本題です。
ここまでで「自由電子近似+ゾンマーフェルト展開」によって化学ポテンシャル(と電子比熱)を計算する準備が整いました。
ここからは、ゾンマーフェルト展開を使わずに化学ポテンシャルを計算することを考えます。

価電子数は温度に依存しないので

n = \int_{-\infty}^{\infty} D(\epsilon)f(\epsilon)\mathrm{d}\epsilon

をμに付いて解く事によって、化学ポテンシャルμをそれぞれの温度で計算することが出来ます。
この式は、解析的には解けないため、Scilabの非線形ソルバfsolveを利用します。

作成したソースコードはChemicalPotential_sce.txtです。

計算結果(自由電子近似+数値計算)


Fig.1に化学ポテンシャルの計算結果を示します。自由電子近似を用いたアルミニウムの場合、ゾンマーフェルト展開を用いた近似と数値積分を利用して求めた結果はほとんど変わらない結果を示します。(高温で若干ずれているのは、本当にゾンマーフェルト展開が妥当でなくなっていっているのか、数値計算の精度の問題かはちょっと自信がないです。計算するエネルギー範囲が狭いせいではなさそう・・・。)

アルミニウムの融点は660.32 ℃(933.47 K)なので、固体の間はゾンマーフェルト展開が破綻することは無いであろう事がわかります。

関連エントリ




参考URL




付録


このエントリで使用したScilab用のソースコードを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 非線形方程式ソルバ 電子比熱 数値積分 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj状態密度モンテカルロ解析常微分方程式トランジスタodeDOSインターフェースPDS5022定電流スイッチング回路分散関係半導体シェルスクリプトレベルシフト乱数HP6632A可変抵抗トランジスタ技術R6452AI2C温度解析ブレッドボードバンドギャップ確率論反強磁性セミナーバンド構造数値積分偏微分方程式非線形方程式ソルバ熱設計絶縁三端子レギュレータISO-I2CA/DコンバータシュミットトリガフォトカプラカオスPWscfGW近似LM358LEDマフィンティン半径発振回路USB数値微分TL431PC817Cサーボアナログスイッチ直流動作点解析補間カレントミラー74HC4053bzqltyチョッパアンプFFT2ちゃんねる開発環境量子力学単振り子電子負荷VESTAQuantumESPRESSO標準ロジックパラメトリック解析ブラべ格子イジング模型アセンブラLDA基本並進ベクトルBSchSMPTLP621失敗談六方最密充填構造コバルト位相図QSGWGGAスイッチト・キャパシタewidth状態方程式VCAキュリー温度繰り返し最適化仮想結晶近似不規則合金熱伝導gfortran相対論抵抗FETMaximaQuantum_ESPRESSOcygwinランダムウォークラプラス方程式スピン軌道相互作用スレーターポーリング曲線マントルシュレディンガー方程式ZnO自動計測QNAP固有値問題ダイヤモンドデータロガー井戸型ポテンシャルTLP552CIFxcrysdenゼーベック係数熱力学条件分岐MCU最小値UPS格子比熱最大値ガイガー管平均場近似過渡解析Writer509スーパーセルFXA-7020ZR差し込みグラフ第一原理計算テスタ起電力OpenMP三角波ubuntuLM555NE555ブラウン運動詰め回路ハーフメタルawkfsolveUbuntuフェルミ面TLP521トランスMAS830LPGACK1026OPA2277フィルタトレーナーバトルEAGLEノコギリ波負帰還安定性ナイキスト線図MBEP-10LMC6622SC1815CapSenseAACircuitPIC16F785入出力固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式ヒストグラム確率論初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型磁気モーメント電荷密度重積分SIC不純物問題擬ポテンシャル状態図cif2cellPWgui二相共存ウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造マテリアルデザインspecx.fフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデル文字列ポケモンGO熱拡散方程式HiLAPW両対数グラフ片対数グラフ陰解法Crank-Nicolson法ifort境界条件連立一次方程式グラフの分割軸ラベルヒストグラム不規則局所モーメントスーパーリーグ円周率Gimp凡例線種シンボルトラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。