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 非線形方程式ソルバ 

comment

Secret

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

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCOPアンプPICCPA強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースDOSPDS5022ecaljスイッチング回路定電流半導体シェルスクリプトレベルシフト乱数HP6632A温度解析ブレッドボードI2CR6452A分散関係トランジスタ技術可変抵抗確率論数値積分反強磁性セミナー非線形方程式ソルバ絶縁バンドギャップ熱設計偏微分方程式バンド構造GW近似カオス三端子レギュレータLEDフォトカプラシュミットトリガISO-I2CA/DコンバータLM358USBカレントミラーTL431マフィンティン半径PC817C数値微分アナログスイッチ発振回路サーボ直流動作点解析74HC40532ちゃんねる標準ロジックチョッパアンプLDAアセンブラFFTbzqltyイジング模型ブラべ格子開発環境補間量子力学電子負荷BSchパラメトリック解析単振り子基本並進ベクトル熱伝導繰り返しGGAMaximaTLP621ewidthSMP相対論抵抗位相図ランダムウォークスピン軌道相互作用六方最密充填構造不規則合金FETコバルト失敗談QSGWcygwinスレーターポーリング曲線スイッチト・キャパシタラプラス方程式gfortranキュリー温度状態方程式条件分岐格子比熱TLP552LM555TLP521三角波NE555過渡解析FXA-7020ZRWriter509テスタ詰め回路MCUマントルダイヤモンドQNAPデータロガーガイガー管自動計測UPS井戸型ポテンシャルawk第一原理計算仮想結晶近似ブラウン運動差し込みグラフ平均場近似fsolve起電力熱力学OpenMPスーパーセル固有値問題最適化最小値VCAシュレディンガー方程式VESTAubuntu最大値面心立方構造PGAOPA2277L10構造非線型方程式ソルバ2SC1815fccフェルミ面等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン正規分布結晶磁気異方性interp1フィルタ初期値ウィグナーザイツ胞c/aルチル構造岩塩構造スワップ領域リジッドバンド模型edeltBaOウルツ鉱構造重積分SIC二相共存ZnOquantumESPRESSOCapSensegnuplotmultiplot全エネルギー固定スピンモーメントFSM合金ノコギリ波フォノンデバイ模型ハーフメタル半金属TeXifortTS-110不規則局所モーメントTS-112等価回路モデルパラメータ・モデルヒストグラムExcel円周率GimpトラックボールPC直流解析入出力文字列マンデルブロ集合キーボードフラクタル化学反応三次元Realforce縮退日本語最小二乗法関数フィッティング疎行列シンボル線種ナイキスト線図陰解法負帰還安定性熱拡散方程式EAGLECrank-Nicolson法連立一次方程式P-10クーロン散乱Ubuntu境界条件MBEHiLAPW軸ラベルトランスCK1026MAS830L凡例PIC16F785LMC662AACircuit両対数グラフ片対数グラフグラフの分割specx.f

最新コメント
リンク

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