数式処理ソフト DERIVE（デライブ） de ドライブ
２９．方程式の数値解法
１．方程式（1変数）の数値解法：ニュートン法
「これまでに4次までの代数方程式は、厳密に解が得られることが分ったのじゃが、5次以上の代数方程式や三角方程式、対数方程式などなどでは、たまたま、因数分解ができるものしか、厳密解が得られない。
そこで、数値的な解法に頼ることになるのじゃが、ここで、グラフは、非常に大切じゃな」
「でもグラフだと、あんまり正確な値が求まらないのじゃないの？」
「もちろん、グラフだけに頼るというわけでもないのじゃが、有力な方法である「ニュートン法」は、グラフで説明すると分かりやすいじゃろう」
「ニュートンって、あの有名なニュートンなの？」
「そうじゃ。
たとえば、x^2－2＝0を計算することで説明しようかの。
答えは、x＝±√2≒±1.4141356・・。
グラフy＝x^2－2を描いてみると、グラフは、y軸に対象で、方程式の正の根は、1と2の間にあることが分る。
ここで、第ゼロ近似値として、x0＝2としてみる。
これから、一般にxn+1＝xn－f(xn）／f'（xn）として計算するのが、ニュートン法じゃ」
「じゃ、とりあえず、やってみるよ。
数式処理ソフト DERIVE（デライブ）で、関数f(n)を再帰的に定義する。
f(n) ：＝ IF(n = 0, 2, f(n － 1)0.5 + 1/f(n － 1))
これで、n=0から10まで計算してみると、こんなふうになったわ。
[2, 1.5, 1.416666666, 1.414215686, 1.414213562, 1.414213562, 1.414213562,
1.414213562, 1.414213562, 1.414213562, 1.414213562]
すごい！ n＝5で6桁ぐらい合っている」
「まったくじゃのう。DERIVEの10桁の標準の精度では、違いが分らない。20桁の精度にしてみると、
[2, 1.5, 1.4166666666666666666, 1.4142156862745098039, 1.4142135623746899106,
1.4142135623730950488, 1.4142135623730950488, 1.4142135623730950488, 1.4142135623730950488,
1.4142135623730950488, 1.4142135623730950488]
真値は、1.4142135623730950488じゃから、n=5で20桁合っていることになるの」
「いったい、どんな原理なのかしら？」
「これが、意外に簡単なのじゃ。関数ｙ＝f(x)として、近似値をx1とすると、x1での接線は、y＝f'(x)x＋bとなる。ここで、接線は（x1，f(x1)）を通るので、
y＝f'(x1)×x＋f(x1)－f'(x1)×x1。この接線がx軸と交わるのは、yがゼロとなることから、次の近似値が得られる。
整理すると、先ほどの逐次近似式が得られるというわけなのじゃ」
「最初の近似値は、粗くてもよいわけね」
「うん。まあ、そうなのじゃ。ニュートン法は、第ゼロ近似値の取り方によっては、発散してしまうこともあるが、収束するときは、普通は、速い」
２．DERIVEでのニュートン法の使い方
「でも、いつも、上のように関数を定義していくのも面倒だわね」
「そう言うと思ったの。ニュートン法は、よく使われるので、次のような便利な方法が用意されているのじゃ。
まずは、「ファイル」メニューの「読み込み」から「Mathファイル」を選択する。そして、EquationSolving.mthというファイルを読み込むのじゃ」
「NEWTON(u, x, x0, n)：＝ ITERATES(x － u/∂(u, x), x, x0, n)
という関数が定義されているのね」
「そうじゃ。さきほどのx^2－2=0という方程式で使い方を説明しよう。
NEWTON(x^2 － 2, x, 2, 10)のように数式入力行に入れるのじゃ。
つまり、uは、解きたい式=0の左辺の関数、xは、解きたい変数を指定する。また、x0は、初期値、nは、繰り返し回数なのじゃな」
「なるほどね。でも、そのMathファイルを読み込まなくても使えるわ」
「あ、本当だな。こりゃ便利だ」
「しっかりしてくださいよ」
「じゃが、DERIVE内部でどのように定義されているかが分るじゃろうが。ヘルプ、読まんでも分るがな」
「おじさん、急に関西弁にならないでください」
「あと、繰り返し回数を指定しない方法もある。
これは、第4引数であるnを省略するのじゃ。
こうすると、DERIVEのオプションで定義されている精度、通常は、10桁じゃが、この範囲に相異なる数列の差が収まるまで繰り返す。この式の場合で10桁では、初期値を除くと5回じゃな」
「解が複素数の場合は、どうなるのかしら？」
「うん。よいところに気がついたの。x^2+1＝0でやってごらん」
「変数を複素数として宣言するのね。でも、初期値を2とすると収束しないわね」
「初期値をたとえば、2+iのような複素数にしてごらん」
「あ、これだと、すぐにiになるわね。でも、当てずっぽうというのは、気持ちが悪いわね」
「そうじゃのう。原理的なことだが、f(z)＝0、ここでzは、複素数としたとき、この近似値は、abs(f(z))を、z＝x+yiとおいて、3次元プロットするか、あるいは、xとyについて、微分するなりして、最小値に近いx、yの組を見つけることじゃが、3次元プロットは、わかりにくいかも知らんな」
「さっきの例だと、元の式が簡単に解けてしまうので、いい例にはならないわね」
「そうじゃな。DOS時代のDERIVEのテキストにこんな例が載っている。
exp（z）－z^2＝0を解け、ここで、zは、複素数じゃ」
「これも、1+i と初期値を置くと、ｚ＝1.588047264 + 1.540223501 iが解として求められる」
３．ITERATE関数とITERATES関数
「ニュートン法の定義でITERATESというのは、何かしら？」
「英語のIterate（繰り返す）から名前をとった関数じゃよ。ま、関数というより手続きといった方が似合っているかも知れん。
ちょっと、まぎらわしいが、ITERATE関数とITERATES関数の2つがあるのじゃよ。
上のITERATES関数は、数列として求める値を関数であり、sの付かないITERATE関数は、n番目の値のみを表示する関数じゃな。
やはり、DOS時代のテキストにこんな例題が載っていた。n番目の素数を求める関数f(n)を定義しなさい。ここで、前に出てきた、NEXT_PRIME（p）を使え」
「f(n)：＝ ITERATE(NEXT_PRIME(k), k, 2, n)でいいわけね」
「惜しいのう。f(n)：＝ ITERATE(NEXT_PRIME(k), k, 2, n－1)じゃな。たとえば、n＝1000は、いくつかな」
「7919だわ」
「うん。NTH_PRIME（n）は、DERIVEの組み込み関数じゃが、同じ値を与えることが分る。
ITERATE関数またはITERATES関数に与える式は、ベクトルでもよいのじゃ。
何回も引用して気が引けるが、同じテキストにフィボナッチ数列の作成例が載っておる。f(0)=0、f(1)=1、f(n)=f(n－1)+f(n－2)じゃ」
「前に再帰関数でやったわね」
「そうじゃ。あのときは、IF関数で再帰的に関数を定義したのじゃが、nが大きくなると時間がかかった。
f(n)：＝ ITERATE([v↓2, v↓1 + v↓2], v, [0, 1], n)とするのじゃが、分るかな？」
「下向き矢印は、ベクトルの要素を表すのだったわね。
ベクトルの2番目の要素を更新していくということなのね。でも、関数の中にベクトルの記号の角括弧が入っているのがわかりにくい」
「確かにのう」
