\documentclass[12pt]{jarticle}
\usepackage{graphicx}
\textwidth=16cm
\oddsidemargin=0cm

\newcommand \Br [1] {\left( #1 \right)}

\renewcommand  \[  {\begin{eqnarray}}
\renewcommand  \]  {\end{eqnarray}}

\begin{document}


\title{工学基礎実験実習 レポート文書作成技術1 \\ | ニュートン法に関する実験 |}
\date{\today}
\author{右田剛史}
\maketitle


\section{はじめに}

方程式の求解は，科学技術計算において頻繁に現れる．解の公式が存在する方
程式では，有限回の四則演算と初等関数によって解を計算できるが，一般の方
程式の解を解析的に表す公式は存在しない．このような場合は，解の近似値を
計算する{\em 数値解法}が有用である．

ニュートン法\cite{numeric,mathwld}は代表的な数値解法である．多くの場
合，ニュートン法は2次の収束を示すため，効率良く解を得ることができる．た
だし，ニュートン法は初期値として解の粗い近似値を与える必要があり，初期
値が適切でない場合，収束までに多くの反復を要したり，収束しない場合もあ
る．また，計算しようとしている解が重解である場合は，収束が遅くな
る．ニュートン法の使用に際しては，このような欠点に留意する必要がある．

以下では， ニュートン法の挙動を理論的に解析し，$f(x)=(x-1)(x+1)^2=0$を
対象として，ニュートン法の性質を調べる実験を行う．


\section{ニュートン法の原理}
\label{sec:bas}

方程式$f(x)=0$の{\em 解}とは，関数$f(x^*)=0$を満たす$x^*$のことを言う．
図\ref{fig:fig1}では，曲線$y=f(x)$と$x$軸が交わっており，この交点の$x$座標が$x^*$である．

いま，解$x^*$の近似値$x_k$が与えられているとする．
点$(x_k,f(x_k))$における曲線$y=f(x)$の接線（図中の斜めの線）
の方程式は$y=f(x_k)+(x-x_k)\cdot f'(x_k)$である．ここで，$f'(x)$は$f(x)$の導関数である．
この接線と$x$軸の交点($x_{k+1}$)は次の式で表される\cite{clo05,mathwld}．
\[
x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}
\label{eq:newton-iter}
\]
図\ref{fig:fig1}の場合，$x_{k+1}$が$x_k$よりも解$x^*$に近い．このため，
適当な初期値$x_0$を与え，式(\ref{eq:newton-iter})によって数列$(x_k)$を
定義すると，この数列は$k\rightarrow\infty$で$x^*$に{\em 収束}することが
期待される．

\begin{figure}[t]
 \begin{minipage}{7.95cm}
  \center
  \includegraphics[scale=.6]{newton-abst.eps}
  \caption{ニュートン法の幾何学的解釈}
  \label{fig:fig1}
 \end{minipage}
 \begin{minipage}{7.95cm}
  \center
  \small
\begin{verbatim}
$ bc -l
x=1.1
x=x-(x-1)*(x+1)*(x+1)/(3*x*x+2*x-1);x
1.00869565217391304348
x=x-(x-1)*(x+1)*(x+1)/(3*x*x+2*x-1);x
1.00007464079119238665
x=x-(x-1)*(x+1)*(x+1)/(3*x*x+2*x-1);x
1.00000000557062401616
x=x-(x-1)*(x+1)*(x+1)/(3*x*x+2*x-1);x
1.00000000000000003104
x=x-(x-1)*(x+1)*(x+1)/(3*x*x+2*x-1);x
1.00000000000000000001
\end{verbatim}
  \caption{実験結果1}
  \label{fig:fig2}
 \end{minipage}
\end{figure}

収束性の議論を厳密にするために，極限値$x^*$との誤差$r_k$を次のように定
義する．
\[
\label{eq:deferr}
r_k := x_k - x^* 
\]
式(\ref{eq:newton-iter})の関数$f(x)$を$x^*$の近傍でテイラー展開して，整
理すると次式が得られる\cite{numeric}．
\[
\nonumber
 r_{k+1} &=& r_k - \left. \frac{r_k f'+ r_k^2 f''/2 + \cdots}{f'+ r_k f'' + \cdots} \right|_{x=x^*}
\label{eq:quadra}
\\      &=& \left. \Br{r_k^2/2} \Br{f''/f'} \right|_{x=x^*} + O(r_k^3)
\]
ただし，$f'(x_k)\ne0$と仮定した．上の式は$r_{k+1}$が$r_k^2$に比例するこ
とを示している．これを2次収束という．これは，$r_k$が十分に小さければ，
式(\ref{eq:newton-iter})の適用によって，正しい桁数がほぼ2倍になることを
示している．逆に，$r_k$が大きいときや$f'(x_k)$が0に近いときには発散する
場合がある．また，$f'(x_k)=0$のときは$x_{k+1}$が計算できない．

なお，重解の場合($f'(x^*)=0$の場合)は 1 次収束である．
また，解付近で2次導関数が 0 になる場合には 3次の収束を示す．


\section{実験}
\label{sec:exp}

ここでは，次の方程式(図\ref{fig:func})にニュートン法を適用して，挙動を
調べる．
\[
f(x)=(x-1)(x+1)^2
\]
解は，$1,-1,-1$である．
導関数は$3x^2+2x-1$であるから，ニュートン法の反復式は次のようになる．
\[
\label{eq:iter}
x_{k+1}=x_k-(x_k-1)(x_k+1)^2/(3x_k^2+2x_k-1)
\]

\begin{figure}[t]
 \begin{minipage}{7.95cm}
  \center
  \includegraphics[scale=.44]{cubic.eps}
  \caption{$(x-1)(x+1)^2$のグラフ}
  \label{fig:func}
 \end{minipage}
 \begin{minipage}{7.95cm}
  \center
  \includegraphics[scale=.44]{ntn.eps}
  \caption{反復回数と誤差の関係}
  \label{fig:fig4}
 \end{minipage}
\end{figure}

初期値を1.1とし，\verb+bc+コマンドを用いてニュートン法の反復を行った結
果を図\ref{fig:fig2}に示す．5回の反復で約20桁まで正しく得られていること
がわかる．初期値が10の場合には11回の反復を要した．

初期値を$-2$とすると重解の$-1$に近づく．しかし，10反復後の$x_{10}$が
$-1.0014\cdots$であり，2桁程度しか正しくない．20桁程度の精度に達するに
は数十解の反復を要すると考えられる．しかし，実際には30数回目以降は桁落
ちのため正しく計算されなかった．

表\ref{tab:表1}に，様々な初期値からニュートン法の反復を10回行った後の値
を示す．$f'(x)=0$となるのは$x=-1,1/3$であるため．$x_k$が$-1$または
$1/3$となってはならない．概ね，初期値が$1/3$を越える場合は1に収束すると
考えられる．また，この場合は2次収束なので，10回以内の反復で収束してい
る．一方，初期値が$1/3$より小さい場合には，重解$-1$に向かうが，1次収束
であるため10回の反復では2桁程度の精度しか得られない．初期値が0の場合は
次に$x_1=-1$となり，それ以上計算できないため N/A (not available) と表示
している．

図\ref{fig:fig4}に，様々な初期値に対する誤差$r_k$の絶対値の推移を対数目
盛で示す．各折れ線付近の数字は初期値を表す． bcの性質上，最後の値は常に
$10^{-20}$の誤差を含むため，折れ線の一番右側が正しくない形に曲がってい
るが，初期値が0.4以上の場合(解1に2次収束する場合)は，ほぼ同じ形状で急速
に誤差が減少していることがわかる．一方，初期値が0.2の場合(解$-1$に1次収
束する場合)は，直線を描いており，この直線が$10^{-20}$に達するには非常に
時間が掛かることが容易に予想できる．

上記の初期値では誤差は単調に減少しているが，一般には誤差が増加すること
もあり得る． $x_k$が$1/3$付近または$-1$付近の値を値をとる場合に
は，$f'(x_k)$が0に近いため$x_{k+1}$が非常に大きな値となり，誤差が増大す
る．このため，ニュートン法の初期値は$x_k$がこのような際どい領域を通過し
ないように選ぶ必要がある．

\begin{table}
 \caption{様々な初期値からのニュートン法の反復結果}
 \label{tab:表1}
 \center
\begin{tabular}{c|c}
\hline
初期値 & 10回の反復後の値 \\
\hline
$-0.2$	& $-0.99964988198316793008$ \\
0	& N/A \\
0.2	& $-1.00357177068625731955$ \\
0.4	& 1.00000000000000000001 \\
0.6	& 1.00000000000000000001 \\
0.8	& 1.00000000000000000001 \\
\hline
 \end{tabular}
\end{table}


\section{まとめ}

ここでは，簡単な方程式にニュートン法を適用し，解と極限値の関係を調べ
た．また，挙動の理論的解析を行った．

実際の応用では，解のわからない難しい方程式を解く必要がある．このような
問題については，今後の更なる検討を要する．



\thebibliography{99}
 \bibitem{numeric} 著者1, 著者2，数値計算の基礎，某出版社，2005．
 \bibitem{clo05} Cox D.A., Little J. and O'Shea D., {\em Using Algebraic Geometry}, Springer, 2005.
 \bibitem{mathwld} \verb|http://mathworld.wolfram.com/NewtonsMethod.html|

\end{document}