[Scilab]関数の宣言とodeの使い方について

お世話になります。

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'); ← φ のプロット

No.1 回答者： cametan_42 回答日時： 2016/01/31 19:12

うーん、これは困りましたねぇ。

誰か詳しい人が回答してくれればイイんですが･･････。

これはまずは質問が錯綜してて、

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となります。)

この手の非線形微分方程式の数値解法だとティピカルな方法論では「４次ルンゲ･クッタ法」なるものが使われるそうです。

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)];
このままでよいとのことでした。

二重振り子のような複雑な例は試していませんが、
独立した振り子を２つ用意した形でシミュレーションしてみたところ、
同じグラフを得られたので恐らく間違っていません。

MATLABの方も興味が出てきたので、紹介してくださった数値解放を含めて勉強してみようと思います。

お礼日時：2016/01/31 20:57