お世話になります。
Scilabの使い方を学んでいる最中で、
以下の振り子を使った数値シミュレーションを参考にしました
「SCILABによる単振り子アニメーション」
http://cacsd2.sakura.ne.jp/2013/11/09/scilab%E3% …
非線形のプログラムの部分で、
function dx=f1(t,x)
th=x(1);
dth=x(2);
dx(1)=dth;
dx(2)=・・・
endfunction
・
・(略)
・
x=[-30/180*%pi; 0]
for i=1:nt, x=[x ode(x:,i),t(i),t(i+1),f1)]; end
・
・(略)
・
x=[x ode(x:,i),t(i),t(i+1),f1)]
の書き方、
・for文で回す必要があるのか
・これは行列の代入を行っているのか
個人的にはここの書き方が分かっていません。
分かりやすく記号書きますと、
恐らく、最初の関数の宣言の部分で、
[θ dθ]という配列を用意して、それに対応するように θ = [θ ode(θ:, i),t(i),t(i+1),f1] という代入をやっているのではないかと思いました。
θだけではなく[θ dθ φ dφ] のようにパラメータが増えたとき、
odeを使ってどのように計算をプログラムで書けばよいか分かりません。
以下のように試しては見ましたが、エラーで怒られてしまい悩んでいます。
function dx=f1(t,x)
th=x(1);
dth=x(2);
phi=x(3);
dphi=x(4);
dx(1)=dth;
dx(2)=・・・
dx(3)=dphi;
dx(4)=・・・
endfunction
・
・(略)
・
x=[-30/180*%pi; 0; 0; 0] ← パラメータが増えた分、初期条件の列も増やす。
for i=1:nt,
x=[x ode(x:,i),t(i),t(i+1),f1) x ode(x:,i),t(i),t(i+1),f1)]; ←ここが分かりません。
end
・
・(略)
・
//プロット
plot(t,180/%pi*x(1,:),'r'); ← θ のプロット
plot(t,180/%pi*x(3,:),'b'); ← φ のプロット
A 回答 (1件)
- 最新から表示
- 回答順に表示
No.1
- 回答日時:
うーん、これは困りましたねぇ。
誰か詳しい人が回答してくれればイイんですが・・・・・・。これはまずは質問が錯綜してて、
1. MATLAB系のプログラミング言語処理系(Octave/SciLab)の使い方に付いて
2. 力学に付いて
3. 非線形微分方程式の数値解法に付いて
の3つが混在してるんですよね。
一つ目。最近はSciLabも人気出てきましたが、多分Octaveのユーザーの方がまだ多いんじゃないか、って事がありますし、また、専門家だと依然として商用であるMATLAB使ってる人が多い事もあって、多分食いつきが悪いと思われます(笑)。
要するに
「俺はMATLABは知ってるけど、SciLabは良く分かんねぇよなぁ」
って事ですね(笑)。
これはGNU Octaveユーザーにしても同じでしょう。
「似てるこたぁ似てる」んだけど、ちと積極的にはなれない(笑)。
2つ目。これは力学の問題であって、実はプログラミング言語処理系のSciLabの問題とも言えない。
高校物理で出てきますが、単振り子の場合、つまり質量mの物体を長さlの紐(棒でもいいですけど)に括りつけて「振り回す」場合、振り幅を角度θとして見ると、運動方程式は
mθ" = -mg*sin(θ)/l
として表現されます。
θが微小角度の場合sin(θ)≒θと「近似」されて
mθ" = -mgθ/l
と表現して構わない。これが参照ページの「線形シミュレーション」ですよね。別にコンピュータ使わなくても解析的に解ける例です。
問題は、「角度が微小じゃない場合」です。これは最初の
mθ" = -mg*sin(θ)/l
をマトモに解かなければならない事を意味しますが、これが解析的には解けないわけですよね。数値的に追っかけてかないといけない。これが件のページでの「非線形微分方程式」になってるわけです。が。
3つ目。件のページ見ても良く分からんです(笑)。いや、だから質問者さんが
「・for文で回す必要があるのか
・これは行列の代入を行っているのか
個人的にはここの書き方が分かっていません。」
と言うのはその通りだよな、です(笑)。
恐らく微小区間内に於いて「常微分方程式のソルバー(つまり問題となってるode)」を適用してる、って意味なんでしょうが・・・。
調べてみる限り、「そんな解法は存在しない」んですよ(笑)。まあ、ザーッと調べた限りの話なんで、ひょっとしたらあるかもしれないですが、少なくとも僕が調べた限りでは「そんな非線形微分方程式の数値解法」はありません。
要するに、一番確実なのは
1. このページの筆者に意図を尋ねる
ないしは
2. 標準的な数値解法のプログラムを自分で書く
しかない、って事です。
(余談ですが、件のページのプログラムは他にも意味が分からないところがあって、例えば、
T=2*%pi*1/sqrt(g/ell)
として大文字Tで変数を定義しています。ちょっと見づらいですが、これは要するに
T=2π/√(g/ell)
なるものを定義してるんですが、実は「このコード内では一回も使われて」いません。使われない変数が何故定義されてるのか意味が分からないです。
また、その前にell=1として定義されてるんですが、だったら
T=2π/√g
でイイんじゃないの?と多分に無意味な事が書かれています。
なお、SciLabはケースセンシティヴ(大文字小文字を判別する)な処理系なんで、当然T≠tとなります。)
この手の非線形微分方程式の数値解法だとティピカルな方法論では「4次ルンゲ・クッタ法」なるものが使われるそうです。
http://slpr.sakura.ne.jp/qp/runge-kutta-ex/
恐らく、これを読んで「自分でSciLabで実装してみる」のが一番確実でしょう。
なお、SciLabのode関数でもルンゲクッタ法が使えるっぽいんで(未確認)、その機能を使ってコードを書いてみても良いかもしれません。
https://help.scilab.org/docs/5.3.0/ja_JP/ode.html
(どうやら、ode関数の第一引数に"rk"と言う文字列を与えると、ルンゲクッタ法で数値解法を試みてくれる模様です。試してみて下さい。)
長く回答をいただきありがとうございます。
>1. このページの筆者に意図を尋ねる
先日、直接電話で問い合わせをしてみたところ、
関数宣言のところでは
th=x(1);
dth=x(2);
phi=x(3);
dphi=x(4);
・
・
・
初期値の部分は
x=[-30/180*%pi; 0; 0; 0] //θ,dθ,φ,dφ
のように単純に増やすだけでよいと、
具体的な説明は聞けませんでしたが、
x=[x ode(x:,i),t(i),t(i+1),f1)];
このままでよいとのことでした。
二重振り子のような複雑な例は試していませんが、
独立した振り子を2つ用意した形でシミュレーションしてみたところ、
同じグラフを得られたので恐らく間違っていません。
MATLABの方も興味が出てきたので、紹介してくださった数値解放を含めて勉強してみようと思います。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- その他(パソコン・スマホ・電化製品) Google ドライブのようにXnBay ストレージ サーバのストレージスペースをコンピュータのエク 2 2023/04/28 19:09
- 船舶・クルーズ Windows10のエクスプローラにて。 1 2022/10/10 20:11
- 電気・ガス・水道 ソーラーパネル初心者です 1 2023/01/01 13:46
- その他(生活家電) アイリスオーヤマの真空パック機 2 2022/12/25 18:26
- 一眼レフカメラ ニコンD7200にシグマ下記レンズはとりつくでしょうか? 1 2023/08/14 18:34
- その他(生活家電) 電熱ヒーターパッド(17×24cm)の電源がすぐに切れるので困っています。 2 2022/12/20 13:31
- サバイバルゲーム このタイプの差し込み口ってなんの種類か分かりますか? 3 2022/07/29 15:32
- ZOZOTOWN このタイプの差し込み口ってなんの種類か分かりますか? 2 2022/07/29 15:31
- その他(自転車) この自転車用ヘルメット、安全なヘルメットではないですか。 10 2023/04/16 07:34
- 一眼レフカメラ ニコンD800はフルサイズですがに下記レンズは付きますでしょうか SIGMA 30mm F1.4 D 2 2023/08/14 19:12
関連するカテゴリからQ&Aを探す
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
VBAのプログラムで、DIAG = 1# ...
-
構造体のデータを丸ごとコピー...
-
ローカル変数の多重定義
-
C言語 構造体の中に共用体を定...
-
VB.NETのStructureというのはど...
-
異なる構造体のデータのコピー
-
Integer変数をカラにしたいので...
-
VBAにてcolorindexを変数に格納...
-
値が代入されてない時
-
typedefをプログラム中で解除す...
-
プログラミング言語の変数と数...
-
文字列の検索&排除をするプロ...
-
C言語のキャストについて
-
構造体の宣言
-
static宣言について
-
変数の初期化について
-
#pragmaについて
-
構造体(問題→回答)
-
Pro*C NUMBER型の...
-
「#undef」と「#define」の使い...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
VBAのプログラムで、DIAG = 1# ...
-
Integer変数をカラにしたいので...
-
C++ 構造体の一括初期化 {0}
-
構造体のデータを丸ごとコピー...
-
C言語 構造体の中に共用体を定...
-
「#undef」と「#define」の使い...
-
VBAにてcolorindexを変数に格納...
-
long型のデータをバイト型の配...
-
値が代入されてない時
-
異なる構造体のデータのコピー
-
構造体のポインタにNULLが入らない
-
VBAの変数のデータ型を変更する...
-
変数の初期化について
-
構造体の初期化方法について
-
ユーザー定義型変数の一括初期化
-
FILE構造体がどのように定...
-
charとucharの違い
-
typedefをプログラム中で解除す...
-
整数から16進数への変換 現在c...
-
VB.NETのStructureというのはど...
おすすめ情報