while文と二分法

1. プログラムの流れの制御(while文)の復習

while文は、継続条件を指定し、不定回数を繰り返しを行います。

構文

while 条件式
    実行コード
end

例: 卒業要件１２４単位と取るまでに、授業を履修して単位と取ること

 total_unit = 0
 while total_unit < 124 
    unit = randi(8);
    total_unit = total_unit + unit;
 end

ここで、randi(n) はnまでの正の整数（即ち、$１,2,\cdots, n$中から疑似乱数を返す。

上記の例では、総単位数が１２４になるまでに、単位を繰り返し取っています。

In [30]:
%以下のコードでは、授業の数を数えている。
 
 total_unit = 0;
 lecture_num=0;
 while total_unit < 124 
    unit = randi(4);
    total_unit = total_unit + unit;
    lecture_num = lecture_num + 1;
 end
 printf("履修科目の数：%d, 履修総単位：%d \n", lecture_num, total_unit)

履修科目の数：49, 履修総単位：127

while文を途中で終了する場合、break文を使用します。

例えば、「履修科目の数は50以上なると卒業できる」という卒業要件がある場合、while文を途中で終了する可能性があります。（実際、このような卒業要件がないです。）

 total_unit = 0;
 lecture_num=0;
 while total_unit < 124 
    unit = randi(4);
    total_unit = total_unit + unit;
    lecture_num = lecture_num + 1;
    if lecture_num >=50
       break
    end
 end

以下のコートを実行してみてください。

In [24]:
total_unit = 0;
lecture_num=0;
while total_unit < 124 
    unit = randi(4);
    total_unit = total_unit + unit;
    lecture_num = lecture_num + 1;
    if lecture_num >=50
       printf("おめでとう！履修科目の数は50以上になるので、卒業でできます。\n");
       break
    end
end
printf("履修科目の数：%d, 履修総単位：%d \n", lecture_num, total_unit)

おめでとう！履修科目の数は50以上になるので、卒業でできます。
履修科目の数：50, 履修総単位：120

演習１

「１からｎまでの整数の２乗の和が50に超える」と満たす最小のｎを求めよ。即ち、 $$ \sum_{i=1}^{n-1} i^2 < 50, ~~ \sum_{i=1}^n i^2 \ge 50 $$ を満たす$n$を計算する。

二分法

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

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

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

注意：

  • 二つのグラフを一つの軸に描くために、「hold on」という命令を使用します。以下の例では、$y=0$と$y=f(x)$を描いています。
In [34]:
%plot -f svg
a=0;
b=3;
x=a:0.001:b;
y=x.^3 - 2;

hold on
grid on
plot(x,y*0,'b-') % x軸を描く
plot(x,y,'r-') % y=f(x)のグラフを描く

p=(a+b)/2
plot(p,p^3-2,'r+') % pにおける関数の値を描く

p =  1.5000
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
  • $a=1,b=3$にすると、$f(a)<0$, $f(b)>0$ が分かります。中間値の定理によって、$a$と$b$の間に$f(x)=0$の解が存在します。
  • $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 [96]:
%plot -f svg
a=0;
b=1.5;
x=a:0.001:b;
y=x.^3 - 2;

hold off
plot(x,y,'r-')
hold on
plot(x,y*0,'b-')

p=(a+b)/2;
plot(p,p^3-2,'ro')
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

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

演習2

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

二分法のアルゴリズム

終了条件:

  • [1)] 許容誤差$TOL$: $(b-a)/2<TOL$（$TOL$は小さい数値とする。例えば、$TOL=1E-5$）
  • [2)] 最大計算回数$N$： $i> N$ ($i$が計算の回数を数える)

計算アルゴリズム

  • 入力：　区間の両端$a$と$b$;許容誤差$TOL$; 最大計算回数$N$。
  • 出力：　近似解$p$或いはエラー情報

詳細

  • [Step 1:] $i=1$, $FA=f(a)$ （初期化）
  • [Step 2:] While ($i\le N$ AND $(b-a)/2 > TOL$) do Step 3-6
  • [Step 3:] Set $p=a+(b-a)/2$, $FP=f(p)$.
  • [Step 4:] If $FP==0$ Then break the while. (Whileループを終了する)
  • [Step 5:] If $FA\cdot FP>0$, Then set $a=p$, $FA=FP$;Else set $b=p$.

  • [Step 6:] Set $i=i+1$.

  • [Step 7:] If $i\le N$ Then Output(p); Else Output("N回の計算で許容誤差の範囲で近似解を見つけなかった。");

演習３：

二分法のアルゴリズムのコードを書いて、$(0,3)$の中で$x^3-2=0$の解を計算しなさい。

ただし、$N=20$,$TOL=0.00001$とする。

In [ ]:
a=0;
b=3;
N=20;
TOL=0.00001;

i=1; FA = a^3-2;
while (...)

end

if i<=N 
  ...
else
  ...
end