床井研究室

大晦日だ

大晦日に私何やってんだろうな。今日はニシン入りそば食えるのかな。この日にあえて何も特別なことをすることなく紅白を見ながら過ごすマンネリとも言える時間を過ごすのは、もしかしたら、こういう当たり前の時間を過ごせること自体がすごく大切でいとおしいことなんだと実感させられた特別な一年だったからかも知れません。またこの一年は、正しいことが一つではないことを思い知らされたことでした。私はいつも「ベストな選択」を求めてきました。というか、「ベストな選択肢」という「正解」が必ずあると信じ込んでいました。でも、そういう「正解」を探しているうちに、自分自身が選択すること、決めるということを先送りにしてきたように思います。「ベストな選択肢」というものは実は存在しなくて、「決めたことにベストを尽くす」ことの方が重要なのではないかと考え直し始めています。ただし、「間違った選択」をしないための注意深さと知恵は身につけておきたいと思います。すみません、ネタが書けませんでした。

年が明けました

あけましておめでとうございます。本年もよろしくお願いします。今年こそはネガティブを返上しようと思います。こういう誰でも見られるところにネガティブなことを書くと、立場上具合の悪いこともある気がします。ちなみに、おみくじは「小吉」でした。「学問：危うし全力を尽くせ」だそうです。まずいです。

微分方程式の数値解法

「ゴムシミュレータ」の力の計算を、ちゃんと物理的な裏付けのあるものにしたいと思います。そのためには、現象を物理的な根拠にもとづいてモデル化する必要があります。しかし、見えている現象をそのまま実験的あるいは経験的にモデル化することが困難な場合もあります。その場合には、現象の非常に短い区間における傾向を記述することにより、うまくモデル化できることがあります。この記述は微分方程式になります。これを解くことにより、現象全体の振る舞いを再現することができます。このあたりの話は私はとっても苦手なので、突っ込まないでください。研究室にはその手の物理の本が何冊か置いてあります。今回は以下のサイトを参考にさせていただきました。ありがとうございます。

• 物理のかぎしっぽ
• LECTURE (Finite Difference Method / Ordinary Differential Equation)
• 計算物理学 (九州大学 野村清英先生)
• 数値解析 (宮城教育大学 瓜生等先生)
• パソコンで数値計算

一階線形常微分方程式

一階線形常微分方程式の基本形は次のようになります。関数 $f$ が現象のある瞬間 $x$ における傾向を表します。

いま、$x$ における $y(x)$ が既知なら、$x + h$ における $y(x + h)$ は、次のようにして求めることができます。

運動方程式

よく知られる $F = ma$ という運動方程式も、運動をモデル化した微分方程式です。空気抵抗を考慮する場合、力を $F(t)$、質量を $m$、加速度を $a$、そして速度を $v$、空気抵抗係数を $k$ とすると、運動方程式は次のようになります。

この力を時刻 $t$ の関数 $F(t)$ とすると、加速度 $a(t, v)$ は次式により求めることができます。

一方、加速度は速度 $v(t)$ を時刻 $t$ について微分して得ることができます。

もし力 $F(t)$ が重力のように時刻 $t$ に関わりなく一定なら、この常微分方程式は解析的に解くことができます。しかし、ここで扱っている「ゴムシミュレータ」では、$F(t)$ が質点の移動に伴って刻々と変化します。$F(t)$ は質点とその近傍との関係で決まるので、これを解析的に解くのは不可能に思えます。そこで、ここではこの微分方程式を数値的に解くことを考えます。現在時刻 $t$ における速度が $v(t)$ であり、時刻 $t$ における加速度が $a(t, v)$ で与えられるとき、$\Delta t$ 後の速度 $v(t + \Delta t)$ は次式により求めることができます。

同様に、位置 $p(t)$ と速度 $v(t)$ との関係は次のようになります。加速度と違って、速度は位置の影響を受けないとします。

したがって、現在時刻 $t$ から $\Delta t$ 後の位置 $p(t + \Delta t)$ は、現在位置を $p(t)$ として次式で求めることができます。

これを図で表すと、次のようになります。

陰解法と陽解法

このような微分方程式の数値解法には、陰解法と陽解法の二つのアプローチがあります。陰解法は区間の後端 $t + \Delta t$ において微分方程式を評価するもので、区間 $\Delta t$ を大きくとっても不安定にならないため、多分こういう用途では主流なんだと思います。布シミュレータでは、Baraff が Large steps in cloth simulation においてこの手法を採用しています。実はうちでも昔こみやーまんがこの手法を使ってました。確かに区間を大きくとっても破綻しないんですが、1 ステップあたりの計算量が結構多くて、こみやーまんは共役勾配法 (CG法) を使ったりしてがんばったものの、インタラクティブ性という点では不満の残る結果でした。

一方、陽解法は微分方程式を区間の前端 $t$ で評価する手法であり、一般には 4 段の陽的 Runge-Kutta 法が使われています。私も遥か昔、寒い中高専の実験室にあった MZ-80K で Runge-Kutta 法の改良手法の Runge-Kutta-Gill 法のプログラムをかじかむ手で作って、結果を放電プリンタにだらだらと印刷して用紙を無駄にしていた思い出があります。

陽解法は区間 $\Delta t$ をむやみに大きく取ると誤差がどんどん累積していって結果が真値から大きく外れて（破綻して）しまいますし、布シミュレータでは Baraff の成果もあるので今更感がないわけでもありません。でも、インタラクティブな応用では 1 ステップの計算量が少ない方が好まれますし、見かけに違和感を感じなければ最終的な誤差は二の次という用途¹では、アニメーションのフレームレートに対する区間長とステップ数を破綻しない程度に調整しておくというアプローチも全く無意味ではないと思うのです。そこで、ここでは陽解法でやってみることにします。

オイラー法

4 段の Runge-Kutta 法は $y(x + h)$ のテイラー展開を 4 次の項までで打ち切って近似するものですが、オイラー法はこれを 1 次の項まで、すなわち $y(t)$ と区間の前端 $t$ における微係数だけを使って $y(x + h)$ を予測（近似）するものです。この方法では、$\Delta t$ 後の速度 $v(t + \Delta t)$ は次式で求めることができます。

$t$ において $v(t)$ は既知ですから、これをそのまま前出の運動方程式に代入すれば、$v(t + \Delta t)$ を求めることができます。同様に位置 $p(t + \Delta t)$ は次式で求めることができます。

この式は $t + \Delta t$ の予測値を $t$ における微係数（傾き）を使って線形予測するものです。それで、この図から見てもわかる通り、この予測は結構外れます。また、外れないように $\Delta t$ を小さくしても、誤差はそれに比例する程度でしか小さくなりません。

改良オイラー法

線形予測ではちょっと外れがでかいので、テイラー展開の 2 次の項までを使って予測（近似）することにします。これには区間の中央の予測値を使う方法と、区間の前端の値と後端の予測値の平均値を使う方法があります。実はどっちがどっちなのかよくわかってないんですけど、多分、前者は改良オイラー法、後者は修正オイラー法と呼ばれているんじゃないかと思います。いずれも 2 段の Runge-Kutta 法に相当するみたいです。

区間の中央値を使う場合、$v(t + \Delta t)$ は次式により求めます。問題になるのは右辺の $v(t + \Delta t / 2)$ の求め方です。これはこの微分方程式そのものの $\Delta t / 2$ 後の解ですから、当然未知の値です。

そこで、この $v(t + \Delta t / 2)$ をオイラー法で予測します。区間が半分になっているので、その分、この予測値の誤差は小さくなります。

改良オイラー法の誤差の大きさは区間の長さの二乗に比例しますから、区間を半分にすれば誤差は四分の一になります。

修正オイラー法

区間の前端の値と後端の予測値の平均値を使う方法があります。

修正オイラー法も誤差の大きさは区間の長さの二乗に比例しますから、区間を半分にすれば誤差は四分の一になります。