3.9.11 発展：関数の入力への過剰依存性

簡単な代数式で与えられる関数は，入力成分が少し変われば出力成分も少し変わる，という性質をよく持つ．しかし，実行手続型の関数には入力成分の変化に過剰に反応するような関数もある．これは，関数で使われる手続が入力の持つ有効数字を累進的に「削っていって」しまうからである．

反復手続に入力値0.1111を与え反復処理を何回か行う．

In[1]:= NestList[FractionalPart[2 #]&, 0.1111, 10]

Out[1]=

値を0.1112に代えて同じ処理を行う．0.1111で見られた結果から累進的に外れていくのが分かる．

In[2]:= NestList[FractionalPart[2 #]&, 0.1112, 10]

Out[2]=

FractionalPart[2 x]の動作は，数 xのバイナリ表記における桁数を考えれば分かりやすい．最初の桁が落とされ，他の桁が左シフトされる．この操作を何回か繰り返すと，より右寄りの桁に対して結果が敏感に反応してもおかしくない．右寄りの桁が xのもとの値にほとんど影響をおよぼさないのである．

FractionalPart[2 x]で行われるシフト操作を再現する． xの最初の8バイナリ桁を見てみる．

In[3]:= RealDigits[Take[%, 5], 2, 8, -1]

Out[3]=

特定の精度で入力すると，事実上桁数が限定されてしまう．1度，桁が「削られ」てしまうと正確な答は求まらない． 精度を上げるためには，もとの入力の数字がもっと分かっていなけらばならないからである． 任意精度の数を使っていれば，Mathematicaは自動的にこのような精度の劣化を見張り， 3.1.5で説明してあるように，によって有効数字が残らないようにした数を示している．

反復処理を繰り返すたびに，累進的に精度が悪化し，ついには桁精度ゼロに陥る．

In[4]:= NestList[FractionalPart[40 #]&, N[1/9, 20], 20]

Out[4]=

各値の桁精度を確認してみる．ゼロ精度は有効桁がないことを示す．

In[5]:= Map[Precision, %]

Out[5]=

厳密に計算すると，答は周期系列になる．

In[6]:= NestList[FractionalPart[40 #]&, 1/9, 10]

Out[6]=

どのような形の近似数でも，それを上記の例のように使うといずれ有効数字がなくなってしまう．任意精度の数を使えば，各ステップの精度が表示されるので少なくとも精度がどう悪化したかが確認できる．しかし，機械精度の数を使うと桁精度が管理できないため，求まった答が意味あるものかどうか分からなくなってしまう．

機械精度の数を使うと桁精度の管理ができない．このため，どう悪化したか確認できない．

In[7]:= NestList[FractionalPart[40 #]&, N[1/9], 20]

Out[7]=

FractionalPart[2 x]の操作を繰り返すたびに，直前の値から有効なバイナリ桁が抽出される．これらの桁の数字がランダムなら（例えば， の桁のように），答もそれに対応してランダムになる．一方，桁の数字が単純なパターンを持つと（例えば，有理数だが），答もそれに対応して単純な形を持つ．

FractionalPart[3/2 x]のような操作を繰り返すと，始めが単純な数であっても，ランダムな数の系列が求まってしまうことがある．この結果は，入力成分への過剰依存の話とは直接関係ないが，実は筆者が1980年代中頃に発見した極めて一般的な現象である．

簡単な入力値だが，一見ランダム的な系列が生成される．

In[8]:= NestList[FractionalPart[3/2 #]&, 1, 15]

Out[8]=

各値が厳密に求まったなら，今度は，安全に数値近似ができる．

In[9]:= N[%]

Out[9]=

厳密数を使い1000回反復し，最後の5回の値を見てみる．

In[10]:= Take[N[NestList[FractionalPart[3/2 #]&, 1, 1000]], -5]

Out[10]=

機械精度で同じ処理をすると，完全に間違った答が求まってしまう．

In[11]:= Take[NestList[FractionalPart[3/2 #]&, 1., 1000], -5]

Out[11]=

反復処理を行う手続きに基づく関数は入力に対して過剰な依存性を示す．似たような状況が微分方程式の解にもいえる．微分方程式の独立パラメータを変化させると，反復処理におけるあるステップから次のステップに移行するときの状況と同じことが起る．

ダフィング(Duffing)の方程式を解く．初期値は1とする．

In[12]:= NDSolve[{x''[t] + 0.15 x'[t] - x[t] + x[t]^3 == 0.3 Cos[t],
x[0] == -1, x'[0] == 1}, x, {t, 0, 50}]

Out[12]=

求まった解をプロットする．

In[13]:= Plot[Evaluate[x[t] /. %], {t, 0, 50}]

Out[13]=

今度は初期値を1.001にして解いてみる．

In[14]:= NDSolve[{x''[t] + 0.15 x'[t] - x[t] + x[t]^3 == 0.3 Cos[t],
x[0] == -1, x'[0] == 1.001}, x, {t, 0, 50}]

Out[14]=

解が前の値から累進的に外れていくのが分かる．

In[15]:= Plot[Evaluate[x[t] /. %], {t, 0, 50}]

Out[15]=