二分法と方程式の根の計算

二分法

区間$(a,b)$上の関数$f$の根を計算するために、二分法が使用されます。

関数$f=x^3-2$を例にして、区間$(0,3)$の中で$f(x)=0$の解の計算方法を説明します。

まず、$y=f(x)$のグラフを書いておきます。

注意：

  • 二つのグラフを一つの軸に描くために、「hold on」という命令を使用します。以下の例では、$y=0$と$y=f(x)$を描いています。
In [17]:
function value=f(x)
    value = x.^3-2;
end
function draw_f(a,b)
    x=a:(b-a)/100:b;
    y=f(x);
    hold on
    grid on
    plot(x,y*0,'b-') % x軸を描く
    plot(x,y,'r-') % y=f(x)のグラフを描く
    p=(a+b)/2;
    plot(p,f(p),'r+') % pにおける関数の値を描く 
end

Step 1.

関数$f$のグラフを描く。

$a=0,b=3$にすると、$f(a)<0$, $f(b)>0$ が分かります。中間値の定理によって、$a$と$b$の間に$f(x)=0$の解が存在します。

In [24]:
a=0;
b=3;
draw_f(a,b)
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -5 0 5 10 15 20 25 0 0.5 1 1.5 2 2.5 3 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a

Step 2

  • 上記のブラフで、$p=(a+b)/2$の値$f(p)>0$が分かります。
  • $f(a)$と$f(p)$の符号が異なるなので、$a$と$p$の間に$f(x)=0$の解が存在します。この場合、$b$を$p$にすると、新しい$(a,b)$で解の計算を進むことができます。
  • 2分法における重要な判定条件式：

$$ f(a)\cdot f(p)<0 $$

In [25]:
%f(a)とf(p)の符号を比較して、pの値を更新します。
p=(a+b)/2;
if f(a)*f(p) < 0
  b = p;
else
  a = p;
end
[a,b]
draw_f(a,b)

ans =

   0.00000   1.50000
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a

Step 3.

新しい$a=0,b=1.5$について、中点$p$で$f(p)<0$ が分かります。また、$f(b)>0$であるので、以下のことが分かる。

$$ f(a)\cdot f(p)>0 \quad \mbox{すなわち}　\quad f(b)\cdot f(p)<0 $$

中間値の定理によって、$p$と$b$の間に$f(x)=0$の解が存在します。

この場合、$a$を$p$にすると、新しい$(a,b)$で解の計算を進むことができます。

In [26]:
%f(a)とf(p)の符号を比較して、pの値を更新します。
p = (a+b)/2;
if f(a)*f(p) < 0
  b = p;
else
  a = p;
end
[a,b]
draw_f(a,b)

ans =

   0.75000   1.50000
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 0.6 0.8 1 1.2 1.4 1.6 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a

演習1

上記の計算を繰り返してすると、解の存在範囲を小さくすることができます。$a$と$b$が十分近いであれば、$(a+b)/2$が$f(x)=0$の近似解となれます。

  • $|b-a|<0.02$まで、上記のコードを繰り返して計算してください。ただし、毎回の計算では、aを更新するか、bを更新するか、調整するのは必要です。

  • 上記のコードをforループに書き直してください。

注意：

ループを終了するとき、breakを使用すること。すなわち、

for k=1:100
   %ここで繰り返し計算のコードを書く。


  if (b-a)< 0.02
    break;
  end
end
In [28]:
% 演習１のコードをここに書いてください。