Scilabユーザーグループ in Japan

初めに †

このページは，Scilabド素人の数学教師TKと，同じくド素人の生徒Aが，Scilabを使って線型常微分方程式を解くまでの道のりを記したものです。全て「フィクション」ですので，余計な詮索はせぬように。

何をしようとしているのか？ †

　そもそものきっかけは，教師TKが「線型常微分方程式の解析解を馬鹿正直に計算するの要する時間と，近似計算をするのに要する時間と，どっちがどれだけ速い(or 遅い)のかを知りたい！」という疑問を持ったことにある。

で，具体的には何をするのか？ †

　教師TKは考えた。無い脳みそを使って考えた結果，次のようなことをしようと決意した。

1. 線型常微分方程式
   
   の解析解を，行列の指数関数を使って
   
   をそのまま計算することで求めたい。ここで
   
   である。
2. その為にはの固有値と固有ベクトルが必須。ただしこれはScilabで簡単に求めることができる。

-->a = [1, 3; 2, 3]

と打ち込めば，行列を入力できる。これを使って固有値をD(対角行列)，固有ベクトルをVに代入するように求めるためには

-->[V, D] = spec(a)

とする。この場合

D =

- 0.6457513 0

0 4.6457513

V =

- 0.8767397 - 0.6354064

0.4809652 - 0.7721779

となる。これはが

というように対角化できていることを表現している。VにはAの固有ベクトル，Dの対角成分にはAの固有値が並ぶ。

これを検算して確認しよう。Scilabのコマンドは

-->V^(-1) * a * V

とすればよい。さすれば

ans =

- 0.6457513 - 3.886D-16

- 2.220D-16 4.6457513

となる。丸め誤差が入るので，正確な対角化はできていないが，固有値・固有ベクトルも良い精度で求められていることが確認できる。

固有多項式を左辺とする代数方程式として解く場合 †

行列の固有値だけ欲しい場合は，代数方程式の根を求めるroots関数を利用する手もある・・・が，近接根の場合は誤差がでかくなる可能性があるので，お勧めはしかねる。

-->a = [1,3;2,3]　 <-- 行列を入力

a =

1. 3.

2. 3.

-->p = poly(a, 'x') <-- 行列aの固有多項式を得る poly(a, 'x', "roots")と同義

p =

2

- 3 - 4x + x

-->roots(p) <-- p = 0を解く

ans =

- 0.6457513

4.6457513

ちなみに，多項式の係数が明示されている場合は

-->p = poly([1,2,3], 'x', "coeff");

-->p

p =

2

1 + 2x + 3x

として与えること。

-->p = poly([1,2,3], 'x');

としてしまうと

-->p

p =

2 3

- 6 + 11x - 6x + x

となって，を展開したものを与えたことになる。その証拠にこれを解かせてみると

-->roots(p)

ans =

1.

2.

3.

↑ほうらね。