方程式の根の近似計算¶

一変数の方程式$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)

x =  1
x =  1.50000000000000
x =  1.41666666666667
x =  1.41421568627451
ans =  1.41421356237310

上記の反復計算では、近似解の数列が得られます。

例えば、以下の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 | ) } $$

以下のレポート課題を解いてください。

計算結果を以下の表にまとめるのはおすすめ。