一変数の方程式$f(x)=0$の根の計算法を考えます。
ニュートン法によって、方程式の近似解を計算することができます。
ニュートン法では、初期値$x_0$をスタートして、$f(x)=0$の解$x$に近づく数列$\{x_n\}$を作ります。 すなわち、
$$ x_0, x_1, x_2, \cdots \rightarrow x $$
まず初めに、予想される真の解に近いと思われる値(初期値$x_0$)をひとつ取ります。
次に、$x_0$におけるグラフの切線と$x$軸の交差点を考えます。この交差点の値は、予想される真の解により近いものとなるのが一般と思われます。以後、この値に対してそこでグラフの切線を考え、同じ操作を繰り返していきます。
以下の計算手順の概念図では、青い線が関数 $f$ のグラフで、その切線を赤で示しています. $x_n$ よりも$x_{n+1}$ のほうが、 $f(x)=0$ の解 x についてのよりよい近似を与えています。
初期値$x_0$から、以下の反復式によって、真の解に収束しているような数列$\{x_n\}$を計算するのは可能です。
$$ x_{n} = x_{n-1} - \frac{ f(x_{n-1}) }{ f^{(1)}(x_{n-1}) } \quad (n\ge 1) $$ ここで、$f^{(1)}(x_{n-1})$は$x_{n-1}$における関数$f$の微分です。
$\sqrt{2}$を算出するために、$f(x)=x^2-2$の解を求めることを考えます。
$$ x^2-2=0 $$
まず, 区間$[0,2]$上の$y=f(x)$のグラフを描きます。
x=0:0.01:2;
y=x.^2-2;
plot(x,y,'b-')
hold on
% (0,0)と(2,0)の間の線分(x軸として)を描画します。
x1=0;y1=0;
x2=2;y2=0;
plot([x1,x2],[y1,y2],'r-')
grid on
この$f$について、ニュートン法の反復計算式は以下のようになります。 $$ x_{n} = x_{n-1} - \frac{ x_{n-1}^2-2 }{2 x_{n-1}} =\frac{1}{2} (x_{n-1} + \frac{2}{x_{n-1}} ) $$
次に、初期値$x_0=1$をスタートして、数列$x_n$を計算します。
% 反復計算のコード
format long %より多くの桁を表示する
x=1
x = 0.5* (x + 2/x)
x = 0.5* (x + 2/x)
x = 0.5* (x + 2/x)
sqrt(2)
上記の反復計算では、近似解の数列が得られます。
例えば、以下のfor文によって、N回反復精算することができます。
N=5; x0=1; x=[x0]; for k=2:N x(k) = 0.5 * ( x(k-1) + 2/x(k-1)) end
実際の計算では、以下の終了条件も使用されます。
$$ |x_n - x_{n-1}|\le 10^{-12} $$
このとき、while文を使って反復計算するのはおすすめです。
また、反復計算では、x_new, x_old二つの変数のみを使って、計算することが可能です。
%上記のfor文の反復計算のコートをwhile文に書き直してください。
ここにコードを書いて下さい。
x_old = 1;
x_new = 0.5*(x_old + 2/x_old );
while abs(x_new-x_old) > 1E-12
... ???
endfor
N=5;
x0=1;
x=[x0];
for k=2:N
x(k) = 0.5 * ( x(k-1) + 2/x(k-1));
end
error_list = abs(x - sqrt(2));
plot(log(error_list),'-+')
grid on
$p$を$f(x)=0$の解または精度の良い近似解とする。反復回数$n$が十分大きい時,反復計算で得られる数列$\{p_i\}$が以下の式を満たすことは考えられる。
$$ \frac{| p_n - p|}{|p_{n-1} - p |^{\alpha}} \approx \frac{| p_{n-1} - p|}{|p_{n-2} - p |^{\alpha}} \approx \lambda $$
ここで、$\lambda$は収束の定数である。
よって,$n\ge 2$のとき,次の式で収束オーダー $\alpha$ を計算できる。一般的に,$n$が大きい時,$\alpha$の計算結果が一定の値となる。
$$ \alpha \approx \frac{ \log (| p_{n-1} - p | / | p_n - p | )}{ \log (| p_{n-2} - p | / | p_{n-1} - p | ) } $$
以下のレポート課題を解いてください。
http://www.ces-alpha.org/download/gs/?label=F5&conf_id=NAP2019&verify=f8237
計算結果を以下の表にまとめるのはおすすめ。
$i$ | Newton method | Order |
---|---|---|
$0$ | $p_0=0.2$ | - |
$1$ | $p_1=?$ | - |
$2$ | $p_2=?$ | $\alpha=?$ ($p_0$, $p_1$, $p_2$から $\alpha$を計算する) |
$3$ | $p_3=?$ | $\alpha=?$ ($p_1$, $p_2$, $p_3$から $\alpha$を計算する) |