初めに

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

何をしようとしているのか?

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

で,具体的には何をするのか?

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

  1. 線型常微分方程式
    \left\{\begin{array}{l}\frac{d\mathbf{x}}{dt} = A\mathbf{x} + \mathbf{g}(t)\\\mathbf{x}(0) = \mathbf{x}_0\end{array}\right.
    の解析解を,行列の指数関数\exp(A)を使って
    \mathbf{x}(t) = \exp(tA)\mathbf{x}_0 + \int^t_0 \exp\[(t - \tau)A\] \mathbf{g}(\tau)d\tau
    をそのまま計算することで求めたい。ここで
    A\in\mathbb{R}^{m\times m},\ \mathbf{x},\ \mathbf{x}_0,\ \mathbf{g(t)}\in\mathbb{R}^m
    である。
  2. その為にはAの固有値と固有ベクトルが必須。ただしこれはScilabで簡単に求めることができる。

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

と打ち込めば,行列A = \left[\begin{array}{cc}1 & 3\\2 & 3\end{array}\right]を入力できる。これを使って固有値をD(対角行列),固有ベクトルをVに代入するように求めるためには

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

とする。この場合

D =

- 0.6457513 0

0 4.6457513

V =

- 0.8767397 - 0.6354064

0.4809652 - 0.7721779

となる。これはA

V^{-1} A V = D

というように対角化できていることを表現している。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

となって,(x-1)(x-2)(x-3)を展開したものを与えたことになる。その証拠にこれを解かせてみると

-->roots(p)

ans =

1.

2.

3.

↑ほうらね。


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-05-04 (日) 21:59:30 (1684d)