初めに †
このページは,Scilabド素人の数学教師TKと,同じくド素人の生徒Aが,Scilabを使って線型常微分方程式を解くまでの道のりを記したものです。全て「フィクション」ですので,余計な詮索はせぬように。
何をしようとしているのか? †
そもそものきっかけは,教師TKが「線型常微分方程式の解析解を馬鹿正直に計算するの要する時間と,近似計算をするのに要する時間と,どっちがどれだけ速い(or 遅い)のかを知りたい!」という疑問を持ったことにある。
で,具体的には何をするのか? †
教師TKは考えた。無い脳みそを使って考えた結果,次のようなことをしようと決意した。
- 線型常微分方程式
の解析解を,行列の指数関数
を使って
をそのまま計算することで求めたい。ここで
である。
- その為には
の固有値と固有ベクトルが必須。ただしこれは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.
↑ほうらね。