方程式の根の近似計算: Secant法

復習

一変数の方程式$f(x)=0$の根を計算するために、Newton法が利用されています。

初期値$x_0$から、以下の反復式によって、真の解に収束しているような数列$\{x_n\}$を計算することは可能です。

$$ x_{n} = x_{n-1} - \frac{ f(x_{n-1}) }{ f'(x_{n-1}) } \quad (n\ge 1) $$

ここで、$f'(x_{n-1})$は$x_{n-1}$における関数$f$の微分です。

セカント法

Newton法の計算には、$f$の微分が使用されています。$f$の微分の計算を避けるために、微分を以下のように近似計算する方法があります。

$$ f'(x_{n-1}) \approx \frac{ f(x_{n-1}) - f(x_{n-2}) }{ x_{n-1} - x_{n-2} } $$

ここで、$f'(x_{n-1})$を算出するために、$x_{n-1}$と$x_{n-2}$二つの近似解を利用することが特徴です。

上記の$f'(x_{n-1})$の近似計算を利用して、セカント法という反復計算法が得られます。

$$ x_{n} = x_{n-1} - \frac{ f(x_{n-1}) ( x_{n-1} - x_{n-2})}{ f(x_{n-1}) - f(x_{n-2}) } \quad (n\ge 2) $$

この反復計算方では、初期値$x_0$と$x_1$が必要です。

計算例

$[0,1.2]$における$f(x)=tan(x^2)-1$の根を計算してくだい。

$$ \tan{x^2}-1=0 $$

まず, 区間$[0,1.2]$上の$y=f(x)$のグラフを描きます。

In [2]:
x=0:0.01:1.2;
y=tan(x.^2)-1;
plot(x,y,'r-')

hold on

 % (0,0)と(2,0)の間の線分(x軸として)を描画します。
 x1=0;y1=0;
 x2=2;y2=0;
plot([x1,x2],[y1,y2],'r-')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -2 0 2 4 6 8 0 0.5 1 1.5 2 gnuplot_plot_1a gnuplot_plot_2a

Step 1. 準備

まず、計算に使用される数学の関数に応じて「MATLABの関数」に作ります。

復習:MATLABの関数の構成 REF:https://jp.mathworks.com/help/matlab/ref/function.html

function 出力変数 = 関数名 (入力変数)

   %与えられる「入力変数」によって、計算行う。
   %最後に計算の結果を「出力変数」にして、結果が返される。

end
In [2]:
function value = myf(x)
    value = tan(x.^2)-1;
end

Step 2. 反復計算式

この$f$について、セカント法の反復計算式は以下のようになります。 $$ x_{n} = x_{n-1} - \frac{ f(x_{n-1}) ( x_{n-1} - x_{n-2})}{ f(x_{n-1}) - f(x_{n-2}) } \quad (n\ge 2) $$

次に、初期値$x_0=0.5$, $x_1=0.6$をスタートして、数列$x_n$を計算します。

In [4]:
% 反復計算のコード
format long  %より多くの桁を表示する
x0=0.5;
x1=0.6;
x2=x1 - myf(x1)*(x1-x0)/(myf(x1)-myf(x0))
x3=x2 - myf(x2)*(x2-x1)/(myf(x2)-myf(x1))
x4=x3 - myf(x3)*(x3-x2)/(myf(x3)-myf(x2))
x5=x4 - myf(x4)*(x4-x3)/(myf(x4)-myf(x3))
x6=x5 - myf(x5)*(x5-x4)/(myf(x5)-myf(x4))
%....
x2 =  1.11511015667020
x3 =  0.725049011967317
x4 =  0.794294713427778
x5 =  0.918416125687412
x6 =  0.879412604736879

Step 3. 反復計算コードの構成

毎回の反復計算では、近似解の列から三つの値(変数)しか関わらないことがありますので、三つの変数で反復計算のコードを書き直します。

以下、変数 x_n_2, x_n_1, xn (それぞれは $x{n-2}$, $x_{n-1}$,$x_n$に対応する)を利用して、コードを書き直します。

$n=6$までコードを書いてください。

In [ ]:
format long  %より多くの桁を表示する

%ーーーーーーーーーーーーーーーー 初期値の設定ーーーーーーーーーーーーーーーーーーー
n=2;
x_n_2 =0.5;
x_n_1 =0.6;
%セカント法により新しい x_n を計算する。

x_n = x_n_1 - myf(x_n_1)*(x_n_1-x_n_2)/(myf(x_n_1)-myf(x_n_2))

%ーーーーーーーーーーーーーーーー 反復計算 1回目  ーーーーーーーーーーーーーーーーーーー

%変数の交代
x_n_2 = x_n_1;
x_n_1 = x_n;

%セカント法により新しい x_n を計算する。
n=3
x_n = x_n_1 - myf(x_n_1)*(x_n_1-x_n_2)/(myf(x_n_1)-myf(x_n_2))

%ーーーーーーーーーーーーーーーー 反復計算 2 回目  ーーーーーーーーーーーーーーーーーーー
% ...

%ーーーーーーーーーーーーーーーー 反復計算 3 回目  ーーーーーーーーーーーーーーーーーーー
% ...

%ーーーーーーーーーーーーーーーー 反復計算 4 回目  ーーーーーーーーーーーーーーーーーーー
% ...

演習1 

上記の反復計算のコートをwhile文に書き直してください。 また、反復計算の終了条件を以下のことを使ってください。

$$ |x_n - x_{n-1}|\le 10^{-13} $$
In [ ]:
%ここにコードを書いて下さい。

n=2;
x_n_2 =0.5;
x_n_1 =0.6;
x_n = x_n_1 - myf(x_n_1)*(x_n_1-x_n_2)/(myf(x_n_1)-myf(x_n_2));
alpha=2;
while (abs(x_n-x_n_1)>1E-13)

  %ここに変数交代のコードを書いてください。
   
  
  %ここに反復計算のコードを書いてください。

  %各Stepの計算結果を出力します。
  % printf("| $%d$ | $x_2=%15.10f $ | $ \\alpha=%5.2f$ | \n",n, x_n, alpha)

end

考察

初期値を$x_0=0.1$, $x_1=0.2$とするとき、反復計算は上手く行けますか?その結果を確認しなさい。

収束オーダーの計算方法

$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 | ) } $$

演習2

セカント法によって、$[0,1.2]$における$f(x)=\tan(x^2)-1$の根を計算してくだい。

初期値:$x_0=0.5$, $x_1=0.6$ 終了条件:$|x_{n}-x_{n-1}|<10^{-13}$

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

$i$ Secant method Order
$0$ $x_0=0.5$ -
$1$ $x_1=0.6$ -
$2$ $x_2=?$ $\alpha=?$ ($a_0$, $a_1$, $a_2$から $\alpha$を計算する)
$3$ $x_3=?$ $\alpha=?$ ($a_1$, $a_2$, $a_3$から $\alpha$を計算する)

演習3

Newton法を使って、$[0,1.2]$における$f(x)=\tan(x^2)-1$の根を計算してくだい。 Newton法とセカント法の収束の速度を比較しなさい。

In [ ]: