8 より進んだ内容(ニュートン法)

4節では、ニュートン法により、実数解を求める方法を学習した。そのとき、連立方程式や複素数解が求められると述べた。そこで、諸君の将来のために、複素数解や連立方程式を求める方法を示しておく。内容は、実数解を求める方法とほとんど同じであるが、少しばかり概念を拡張する必要がある。いままでの学習内容を十分理解していれば、これから述べることは分かるであろう。理解できないにしても、いずれは理解できるものと期待している。内容が理解できなくても、なにか面白そうだと思ってもらえれば、十分である。さあ、実数解が求められた者は、より難しい非線形方程式の近似解を求めよう。

8.1 非線型方程式の複素数解

8.1.1 実数解の場合の復習

複素数に進む前に実数のニュートン法の復習を行う。4節ではでは、図によるニュートン法の説明を行い、漸化式を示した。この方法は、直感的で分かりやすいが、複素数解や非線形連立方程式を考える場合は無理がある。図で示すことが出来るのは2次元までで、それ以上になると図に示すことが不可能であるからである。そこで、別のアプローチから漸化式を導くことにする。異なった説明をしているようであるが、その根本精神は全く同じで、

解に近いところでは、直線で近似できる⁶。

ということである。

それでは、漸化式を求めることにする。いつものように、の方程式の解を $x=\alpha$ とする。即ち、 $f(\alpha)=0$ である。そして、番目の近似解をとする。ここから、 $\Delta x$ だけ移動したところの値は、

$\begin{equation*}\begin{aligned}f(x_i+\Delta x)&=f(x_i) +f^\prime(x_i)\Delta x}\\ &\simeq f(x_i)+f^\prime(x_i)\Delta x \end{aligned}\end{equation*}$

となる。もし、 $f(x_i+\Delta x)=0$ 、即ち、 $\alpha=x_i+\Delta x$ となるように、 $\Delta x$ を選ぶことができたら、解の計算は簡単である。この場合、式 (15)の最後の式から、

$\displaystyle \Delta x \simeq -\frac{f(x_i)}{f^\prime(x_i)}$

(16)

となる。したがって、 $\alpha=x_i+\Delta x$ から、次の近似解は

$\displaystyle x_{i+1}=x_i-\frac{f(x_i)}{f^\prime(x_i)}$

(17)

となる。前回、図により求めた漸化式と同じである。異なる説明であったが、内容はまったく同じであることを理解して欲しい。

8.1.2 複素数解の場合

8.1.2.1 漸化式

実数とまったく同じ議論が、複素数でも成り立つ。ただし、複素関数で重要な特異点付近では、この方法は使えない。テイラー級数ではなく、ローラン級数の

乗よりも小さい項が重要となるからである。実数の場合も、特異点はだめなのと同じである。

実数とまったく同じ議論より、方程式

$\displaystyle w(z)=0$

(18)

の近似解は、漸化式

$\displaystyle z_{i+1}=z_i-\frac{w(z_i)}{w^\prime(z_i)}$

(19)

より求めることができる。この式の算出は、先ほどの実数の場合と全く同じである。

前回とは異なり、実数の場合の漸化式をグラフを用いないで説明したのは、複素数に拡張するためである。複素数のグラフは大変である。

8.1.2.2 プログラム作成のために

つい最近まで、FORTRANと違ってC言語では複素数をそのまま扱うことができなかった。これがC言語の弱点になっていたが、現在ではそれが克服された。FORTRAN同様に複素数もそのまま扱えるのである。これは非常にありがたい。

C言語で複素数を使うためには、教科書p.471に書かれているようにすればよい。すなわち、

	#include <complex.h>

とヘッダーファイルをインクルードし、変数の型を

	float complex
	double complex
	long double complex

のようにする。通常は、double complexを使うこと。

四則演算は特に気にすることもなく、普通に演算子(+,-,*,/)が使える。また、表1のような関数が用意されている。

**表 1:** `complex.h`で定義されている関数。関数の引数と戻り値は同じ型である。
関数名	倍精度	単精度	拡張倍精度
三角関数	`csin()`	`csinf()`	`csinl()`
	`ccos()`	`ccosf()`	`ccosl()`
	`ctan()`	`ctanf()`	`ctanl()`
逆三角関数	`casin()`	`casinf()`	`casinl()`
	`cacos()`	`cacosf()`	`cacosl()`
	`catan()`	`catanf()`	`catanl()`
双曲線関数	`csinh()`	`csinhf()`	`csinhl()`
	`ccosh()`	`ccoshf()`	`ccoshl()`
	`ctanh()`	`ctanhf()`	`ctanhl()`
逆双曲線関数	`casinh()`	`casinhf()`	`casinhl()`
	`cacosh()`	`cacoshf()`	`cacoshl()`
	`catanh()`	`catanhf()`	`catanhl()`
指数関数	`cexp()`	`cexpf()`	`cexpl()`
自然対数	`clog()`	`clogf()`	`clogl()`
絶対値	`cabs()`	`cabsf()`	`cabsl()`
平方根	`csqrt()`	`csqrtf()`	`csqrtl()`
べき乗	`cpow()`	`cpowf()`	`cpowl()`
実部	`creal()`	`crealf()`	`creall()`
虚部	`cimag()`	`cimagf()`	`cimagl()`
偏角	`carg()`	`cargf()`	`cargl()`
複素共役	`conj()`	`conjf()`	`conjl()`
リーマン球の射影	`cproj()`	`cprojf()`	`cprojl()`

8.2 非線型連立方程式の実数解(2元の場合)

前節では、ニュートン法による複素数の近似解を求める方法を示した。「非線型連立方程式」の近似解が求めれれば、概ねニュートン法の学習は終わりである。ちょっと難しいが、ニュートン法の学習の仕上げとして、「非線型連立方程式」の実数解を求める方法を示す。非線型連立方程式の複素数解を求めることが残っているが、この講義では示さない。今までのことを理解していれば、その方法も直ぐに理解できるであろう。興味のある人、あるいは必要に迫られた人は、自分で計算方法を考えてみよう。

8.2.1 非線型連立方程式とは

今まで、諸君は、「非線型の方程式」あるいは「線形の連立方程式」は解いたことがある。例えば、前者は、

$\displaystyle x^2-3x+2=0$

(20)

のようなものである。後者の例は、

$\begin{equation*}\left\{ \begin{aligned}3x+2y+z&=10\\ x+y+z&=6\\ x+2y+z&=11 \end{aligned} \right.\end{equation*}$

である。非線型方程式は未知数が2次以上のものをいい、連立方程式は未知数が2個以上のものをいうのである。非線型とは、直線でないという意味である。未知数が2次以上のもの、例えば

が式に含まれると、それは直線にならないので、非線型方程式になる。直線でないという意味からも、 $\sin x$ も非線型方程式を形づくる。この場合、

の次数は無限である。

複素数解まで含めると、非線型なn次方程式にはn個の解がある。線形な連立方程式、n元1次方程式の場合、係数行列が特異でない限り、1個の解がある。では、非線形な連立方程式、n元 m次方程式の場合、複素数を含めた解の数はm個のように思えるが、正しいのだろうか?。数学の先生に聞くと、正しいということである。

また、方程式の数と未知数の数は一致しなくてはならないのは、通常の連立方程式と同じである。それらの数が同じでも、線形連立方程式では、係数行列の行列式がゼロの場合、解は一意に決まらない。非線型の連立方程式の場合、これはどのような場合に対応するのだろうか?。私には、分からない。かなり難しく興味深い問題のように思えるが、ここではそのことは考えないことにする。

8.2.2 非線型連立方程式の例

実例を使って、計算方法を示す。例として

$\begin{equation*}\left\{ \begin{aligned}&(x-3)^2+y^2-3=0\\ &\sin x+e^{y-1}-1=0 \end{aligned} \right.\end{equation*}$

の近似解を求めることを考える。2元?次非線型連立方程式である。?次とは、いささかいい加減に書いているが、勘弁してもらいたい。無限次といってよいような気がするが自信が無いので?マークを付けておく。

さて、この方程式の解であるが、それをグラフに示す。2元であればグラフに書くことができるのである。以下の議論は、任意の元の方程式でも成り立つことは理解して欲しい。これらの方程式をグラフに書くと図11のようになる。図中に示すように、点AとBに実数解があるのが分かるであろう。

**図 11:** 非線型方程式のグラフと実数解
$\includegraphics[keepaspectratio, scale=1.0]{figure/nonlinear_eqs/ZeroPoint.eps}$

初期値から出発して、解であるAやB点に近づく方法を考える。そこで、次のような関数を考える。

$\displaystyle f(x,y)$	$\displaystyle =(x-3)^2+y^2-3$	(23)
$\displaystyle g(x,y)$	$\displaystyle =\sin x+e^{y-1}-1$	(24)

もちろん、

と

が同時に成り立つ、

を求めたいわけである。

いつものように、この非線型連立方程式の解を $(\alpha_x, \alpha_y)$ とする。当然、 $f(\alpha_x, \alpha_y)=0$ かつ $g(\alpha_x, \alpha_y)=0$ である。そして、番目の近似解をとする。ここから、 $(\Delta x, \Delta y)$ だけ移動したところの値は、

$\begin{equation*}\begin{aligned}f(x_i+\Delta x, y_i+\Delta y)&=f(x_i, y_i) +\fra... ...l x}\Delta x +\frac{\partial f}{\partial y}\Delta y \end{aligned}\end{equation*}$

となる。

の場合も全く同じである。それら、2つをまとめ、 $\simeq$ を

に直すと

$\displaystyle f(x_i+\Delta x, y_i+\Delta y)$	$\displaystyle =f(x_i, y_i) +\frac{\partial f}{\partial x}\Delta x +\frac{\partial f}{\partial y}\Delta y$	(26)
$\displaystyle g(x_i+\Delta x, y_i+\Delta y)$	$\displaystyle =g(x_i, y_i) +\frac{\partial g}{\partial x}\Delta x +\frac{\partial g}{\partial y}\Delta y$	(27)

となる。

ここで、 $f(x_i+\Delta x, y_i+\Delta y)=0$ かつ $g(x_i+\Delta x, y_i+\Delta y)=0$ となるように、 $\Delta x$ と $\Delta y$ を選ぶとする。このようにするためには、 $\Delta x$ と $\Delta y$ はつぎの連立方程式を満たせばよい。式(26)と (27)の左辺をゼロとおき式を整理すれば

$\displaystyle \begin{pmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{... ...\Delta y \end{pmatrix} = \begin{pmatrix}-f(x_i,y_i)\\ -g(x_i,y_i) \end{pmatrix}$

(28)

となる。この連立方程式を解いて、 $(\Delta x, \Delta y)$ を求める。 $\alpha_x \simeq x_i+\Delta x$ したがって、 $\alpha_y \simeq x_i+\Delta x$ から、次の近似解は

$\begin{equation*}\left\{ \begin{aligned}x_{i+1}&=x_i+\Delta x\\ y_{i+1}&=y_i+\Delta y \end{aligned} \right.\end{equation*}$

となる。これが、非線型連立方程式の漸化式である。

8.2.3 連立で無い場合とのアナロジー

以前の授業で示した方程式の実数解や複素数解を求めたのと同じようなことを、ここでも行った。式も似ているし、考え方も同じである。以前は、

解に近いところでは、直線で近似できる⁷。

のような性質を利用した。2元の非線型連立方程式でも同じで、

解に近いところでは、平面で近似できる⁸。

という性質を利用している。この性質を定量的に表したものがテイラー展開である。なるほどテイラー展開は便利なものである。

8.3 非線型連立方程式の解(多元の場合)

前章では、2元の非線型連立方程式のニュートン法での計算方法を示した。ここでは、それを一般化する。ここで示す方法は、複素数解にも適用できる。

N元の非線型連立方程式は、

$\begin{equation*}\left\{ \begin{aligned}&f_1(x_1+x_2+x_3+\cdots+x_N)=0\\ &f_2(x_... ...0mm}\vdots\\ &f_N(x_1+x_2+x_3+\cdots+x_N)=0 \end{aligned} \right.\end{equation*}$

と書くことができる。未知数は、

$\displaystyle \boldsymbol{X}=(x_1,x_2,x_3,\cdots+x_N)$

(31)

とベクトルで表現する。すると、

番目の方程式は、 $f_i(\boldsymbol{X})$ と書き表されるので、表現が簡単になる。これを、先ほどと同じようにテイラー展開すると

$\displaystyle f_i(\boldsymbol{X}+\Delta\boldsymbol{X})= f_i(\boldsymbol{X})+ \f... ...\cdots+ \frac{\partial f_i}{\partial x_N}\Delta x_N+ O(\Delta \boldsymbol{X}^2)$

(32)

となる。 $i=1,2,3,\cdots,N$ の全てににおいて、 $f_i(\boldsymbol{X}+\Delta\boldsymbol{X})=0$ になるように、 $\Delta\boldsymbol{X}=(\Delta x_1, \Delta x_2, \Delta x_3,\cdots ,\Delta x_N)$ を選ぶ。そのように選ぶためには、2次以降の高次の項を無視すると

$\displaystyle \begin{pmatrix}\frac{\partial f_1}{\partial x_1} & \frac{\partial... ...bol{X}) \\ -f_3(\boldsymbol{X}) \\ \ldots \\ -f_N(\boldsymbol{X}) \end{pmatrix}$

(33)

の線形であるN元1次連立方程式が成り立つ。これを解いて、 $\Delta\boldsymbol{X}=(\Delta x_1, \Delta x_2, \Delta x_3,\cdots ,\Delta x_N)$ を求める。そうすると、より真の解に近い $\boldsymbol{X}^{new}$ は、 $\boldsymbol{X}^{new}=\boldsymbol{X}^{old}+\Delta \boldsymbol{X}$ と計算できる。しつこいようであるが、成分で書き表すと

$\begin{equation*}\left\{ \begin{aligned}x_1^{new}&=x_1^{old}+\Delta x_1 \\ x_2^{... ...\ &\vdots\\ x_N^{new}&=x_N^{old}+\Delta x_N \end{aligned} \right.\end{equation*}$

となる

非線型の連立方程式を線形の連立方程式で計算しているわけである。解きやすい式になった分、反復計算が必要となっている。

ホームページ: Yamamoto's laboratory
著者: 山本昌志

Yamamoto Masashi
平成19年6月24日