123456789

「数理物理」講義ノート

第四章：数値微分と数値積分

コンピューターで微分をしよう

微分の概念そのものは、数学者が考案されたかのように思われがちだが、かの物理学者ニュートンによって考え出されたものである。速度の考え方は、単位時間当たりに進んだ距離であるが、この時単位時間というのを限りなく微小な時間Δｔにもっていくことに微分の本質がある。速度をVとするとこれは、刻々変化していくから時間ｔの関数V(t)と考えることができる。速度は微分を使って、X(t)を距離として、

と表せる。「Δtを限りなく０に近づける」という操作は、いかにも数学的である。これを、元々の物理的な考え方に戻せば、どうなるか。

のようになり、近似的に求まることになる。「十分ちいさな値」は、人間が用意するには大変であるが、コンピューターならできる。では、どのくらい小さな値でも良いのか？一章でも述べたように、コンピューターも有限の賜物であることを忘れてはならない。それは、一般にXの関数の性質による。実際に色々な初等関数をXに見立てて、「微分」してみよう。以下のプログラムでは、例として　(t)=, (t)=, (t)=sin(t) についてｔ＝１で計算を行った。

これを実行すれば、以下のようになる。

dt=0.100000 15.937431 10.000000 2.858844 2.718282 0.497364 0.540302
dt=0.010000 10.462201 10.000000 2.731919 2.718282 0.536090 0.540302
dt=0.001000 10.045647 10.000000 2.719879 2.718282 0.539958 0.540302
dt=0.000100 10.006428 10.000000 2.720356 2.718282 0.540614 0.540302
dt=0.000010 10.013581 10.000000 2.741814 2.718282 0.542402 0.540302
dt=0.000001 9.536743 10.000000 2.622604 2.718282 0.536442 0.540302
dt=0.000000 11.920929 10.000000 4.768372 2.718282 1.192093 0.540302
dt=0.000000 0.000000 10.000000 0.000000 2.718282 0.000000 0.540302
dt=0.000000 0.000000 10.000000 0.000000 2.718282 0.000000 0.540302
dt=0.000000 0.000000 10.000000 0.000000 2.718282 0.000000 0.540302

数字は、と並んでいる。このように、単精度の計算では、「うまく」Δｔを選ばないと良い結果が得られない。

次に、それでは倍精度で計算すればどうか？プログラム中の「float」の部分を全て「double」と書き直して実行すればよい。その結果は、期待していたように、

dt=0.100000 15.937425 10.000000 2.858842 2.718282 0.497364 0.540302
dt=0.010000 10.462213 10.000000 2.731919 2.718282 0.536086 0.540302
dt=0.001000 10.045120 10.000000 2.719641 2.718282 0.539881 0.540302
dt=0.000100 10.004501 10.000000 2.718418 2.718282 0.540260 0.540302
dt=0.000010 10.000450 10.000000 2.718295 2.718282 0.540298 0.540302
dt=0.000001 10.000045 10.000000 2.718283 2.718282 0.540302 0.540302
dt=0.000000 10.000005 10.000000 2.718282 2.718282 0.540302 0.540302
dt=0.000000 10.000000 10.000000 2.718282 2.718282 0.540302 0.540302
dt=0.000000 10.000001 10.000000 2.718282 2.718282 0.540302 0.540302
dt=0.000000 10.000001 10.000000 2.718283 2.718282 0.540302 0.540302

良い結果を得る。しかし、Δｔをいくらでも小さく取れるとは限らないことを忘れないように。答えが正しいかどうかを知る簡単なチェックは、Δｔを少し変化さてみると分かる。

コンピューターで積分をしよう

積分は微分の逆の演算であることを考え、上で行ったことをもとに、式を立てていこう。

これを式変形して、

と書ける。これは、漸化式になっている。一般化すれば、

となり、任意のｔについて求まる。このような方法を線形近似法または、オイラー法という。精度はΔｔを小さくすれば良いが、和を多く行わなくてはならなくなる。X(T)をもとめたいとすれば、

ｎ回行う必要がある。(CPU-TIMEがかかる）V(t)=, V(t)=およびV(t)=sin(t)の場合、X(1)=1の初期値から、X(2)を計算した。C言語での、プログラムは以下のようになる。

結果は、以下のようになる。

dt=0.100000 140.186573 187.090909 5.441127 5.670774 1.952261 1.956449
dt=0.010000 182.018491 187.090909 5.647459 5.670774 1.956102 1.956449
dt=0.001000 186.579835 187.090909 5.668439 5.670774 1.956415 1.956449
dt=0.000100 187.039763 187.090909 5.670541 5.670774 1.956446 1.956449
dt=0.000010 187.085794 187.090909 5.670751 5.670774 1.956449 1.956449
dt=0.000001 187.090398 187.090909 5.670772 5.670774 1.956449 1.956449
dt=0.000000 187.090858 187.090909 5.670774 5.670774 1.956449 1.956449
dt=0.000000 187.090904 187.090909 5.670774 5.670774 1.956449 1.956449
dt=0.000000 187.090909 187.090909 5.670774 5.670774 1.956449 1.956449

このように関数系によって、Δｔの大きさがまちまちである。微分の時と同じく収束しているかのチェックが大切である。
９行目のばあいΔｔ＝０．００００００００１であるから、和の計算のループは、１０億回行ってることになる。積分和をほんの２０回ぐらいで同じ精度の計算をすることができる方法がある。その方法の一つにガウス・ルジャンドル積分法がある。（後述）

この他、等間隔に関数を求め、積分する方法に台形公式やシンプソンの公式がある。

「科学技術は、補間だ！」

関数をモデルに考えてみよう。初等関数（、sin(x) ,exp, log など）で表せる関数ならば、ｘのとりえる範囲で、自由に計算できる。しかし、現実の問題は、大抵、データと呼ばれる離散型の数値で与えられている。
例えば、｛x、f(x)｝のデータを以下のように与えて、プロットさせると、[mathematica]

このように描ける。「補間」とは、この場合、データに載っていない値を予想することである。それには、２種類の補間がある。
•変数の区間内の補間　⇒　内挿補間　　(interpolation)
•変数の区間外の補間　⇒　外挿補間　　(extrapolation)
例の場合、区間とは１<ｘ<１０である。これを、私は「知られている世界（既知領域）」、また、区間外（ｘ<１とｘ>１０）は「知らない世界（未知領域）」と考えたい｡ そのように考えると、内挿補間をより詳しく予想していくことは、「技術」で、外挿していく作業は、「科学」をすることに類似性を感じとれる。パラダイム１　パラダイム２

例の場合を考えると、ｘ＝３．５のｆ（ｘ）の内挿は、ｆ（３）とｆ（４）が知られているから、それは丁度真中の値であることから、ｆ（３．５）~（ｆ（３）＋ｆ（４））＝と「予想」する。しかし、ｘ＝１２の値は、簡単には予想できない｡

皆さんは、既に想像できたと思いますが、これはｆ（ｘ）＝が考えられる。これは「仮説」で、これが正しければ、ｘ＝１２は当然、ｆ（１２）＝１４４と求まる。また、ｆ（３．５）＝と計算ができる。では、上のは、全くの間違いだったのかというと、そうではない。「誤差」を計算すると、

２．０４％になり、そんなに悪くない。なぜ悪かったのかというと、たった２つのデータ（ｆ（３）とｆ（４））しか使わなかったからである。多くのデータ点を使って、精度の良い内挿を行う方法が知られている。そのなかには、スプライン内挿法、ラグランジュ内挿法が代表的である。ラグランジュ内挿法について考えていこう。

ラグランジュ内挿法

この補間法は、関数を多項式で表すことを仮定する。これは自然な仮定で、一般に正則な関数は、テーラー展開できるからである。

つまり、データの点｛｝（i=1,...,N）で、f(x)がデータの値｛｝を満たす、

を要請すれば、係数のかずは、N個で、M=N-1となる。（未知数の数と式の数）ここで、係数を求めることを考えずに、或る多項式を用いて、

を仮定すれば、f(を満たすためには、

であれば良い。（但し、＝１(i=jの場合)、＝０(その他の場合)　）これを満たすN-1次の多項式は、ユニークに与えられる。

N-1次式になるのは、分子がないことに注意。これを、ラグランジュの積関数と呼ぶ。

ガウス・ルジャンドル積分法

定積分の計算をしよう。変数変換すれば、任意の積分区間を -1<t<1とすることができる。（ x⇒　 ）

上のラグランジュの補間の考え方から、関数f(x)を、ラグランジュの積関数とで、

と書けるので、積分 I は、

と書き換わる。ここではを積分の「重み」という。この重みは、積分点が既知ならば、予め求めておくことができる。

シンプソンの公式の証明：

連続する３点を２次関数で近似して、その解析的な積分値を積み重ね加えていく方法。
３点を(−１,)(０,),(１,)とすれば、それぞれの点での、、がもとまる。

問い）４点を選び、３次関数で近似する場合はどうか？

次に、積分を行う前に、積分点｛｝を決めることができる場合を考える。今、N次の関数を仮定すると、

ラグランジュ積関数は、

と書ける。ここで、をルジャンドル関数という特殊関数に選び、積分点｛｝をそのゼロ点（＝０）に選ぶと、重みは解析的に求まる。（厳密には、ルジャンドル関数をとおくには定数倍違うが、ラグランジュ積関数は変わらない。）

問題：ガウス・ルジャンドル積分のプログラム

ルジャンドル関数のゼロ点と重みは、N=２，４，６について、次の様に与えられている。
ｎ　　
2    +-  0.57735027      1.00000000  
4    +-  0.86113631      0.34785458
      +-  0.33998104      0.65214515
6    +-  0.93246951      0.17132449
      +-  0.66120939      0.36076157
      +-  0.23861919      0.46791393      

これらを用いて、オイラー法で行なった、V(t)=, V(t)=およびV(t)=sin(t)について

の計算を求めよ。

第五章へ　　 第三章へ

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　H.Kamada

---

Converted by Mathematica March 15, 2002