$b-a<0.02$ となるまで、二分法の操作を手動で繰り返し行ってください。「手動で行う」演習ですので、if文、for文、while文を使わないこと。Markdownとして説明を付ける必要はなく、コードのみを書けばよいです。
function value = f(x)
value = x.^3-2;
end
a = 0;
b = 3;
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
b = p; %f(a)*f(p)<0なので、aとpの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
a = p; %f(a)*f(p)>0なので、pとbの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
a = p; %f(a)*f(p)>0なので、pとbの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
b = p; %f(a)*f(p)<0なので、aとpの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
a = p; %f(a)*f(p)>0なので、pとbの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
b = p; %f(a)*f(p)<0なので、aとpの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
a = p; %f(a)*f(p)>0なので、pとbの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
a = p; %f(a)*f(p)>0なので、pとbの間に厳密解が存在
p = (a+b)/2;
printf("a=%f, b=%f, b-a=%f, p=%f, f(p)=%f, f(a)*f(p)=%f\n",a,b,b-a,p,f(p),f(a)*f(p))
%b-a<0.02となったので、終了
a=0.000000, b=3.000000, b-a=3.000000, p=1.500000, f(p)=1.375000, f(a)*f(p)=-2.750000 a=0.000000, b=1.500000, b-a=1.500000, p=0.750000, f(p)=-1.578125, f(a)*f(p)=3.156250 a=0.750000, b=1.500000, b-a=0.750000, p=1.125000, f(p)=-0.576172, f(a)*f(p)=0.909271 a=1.125000, b=1.500000, b-a=0.375000, p=1.312500, f(p)=0.260986, f(a)*f(p)=-0.150373 a=1.125000, b=1.312500, b-a=0.187500, p=1.218750, f(p)=-0.189728, f(a)*f(p)=0.109316 a=1.218750, b=1.312500, b-a=0.093750, p=1.265625, f(p)=0.027287, f(a)*f(p)=-0.005177 a=1.218750, b=1.265625, b-a=0.046875, p=1.242188, f(p)=-0.083268, f(a)*f(p)=0.015798 a=1.242188, b=1.265625, b-a=0.023438, p=1.253906, f(p)=-0.028507, f(a)*f(p)=0.002374 a=1.253906, b=1.265625, b-a=0.011719, p=1.259766, f(p)=-0.000740, f(a)*f(p)=0.000021
二分法のアルゴリズムを実現するコードを書き、関数 $f(x)=x^3-2$ に対して方程式 $f(x)=0$ の区間 $[0,3]$における近似解を求めてください。ただし、許容誤差 $TOL=0.00001$ とします。
function value = f(x)
value = x.^3-2;
end
a = 0;
b = 3;
TOL = 0.00001;
n = 0;
while b-a >= 2*TOL
n = n+1;
p = (a+b)/2;
if f(a)*f(p) < 0
b = p;
else
a = p;
end
end
p = (a+b)/2;
printf("反復回数n=%d, 近似解p=%f\n",n,p)
反復回数n=18, 近似解p=1.259920
求めた近似解を厳密解と比較してみると、誤差は確かに許容誤差 $TOL$ よりも小さくなっていることが分かります。
q = 2^(1/3); %厳密解
printf("近似解=%.10f, 厳密解=%.10f, 近似解と厳密解の誤差=%.10f\n",p,q,abs(p-q)) %printfでは小数点以下を指定桁数で表示することが可能
近似解=1.2599201202, 厳密解=1.2599210499, 近似解と厳密解の誤差=0.0000009297
二分法のアルゴリズムにおける区間の両端の初期値 $a$ と $b$、許容誤差 $TOL$、反復回数 $n$ の関係について考察してください。Markdownのセルに解答を書いてください。
二分法のアルゴリズムにおいては一回の反復で区間の幅が半分になるので、$n$ 回反復した後の区間の幅は
$$ \frac{b-a}{2^n} $$です。反復の終了条件を考えると、
$$ \frac{b-a}{2^n}<2\cdot TOL $$が成り立つことが分かります。これを $n$ について解くと、
$$ n>\log_2 \frac{b-a}{TOL}-1 $$が得られます。つまり、 反復回数 $n$ はこの不等式の右辺より大きな最小の自然数になります。
実際に $a=0$、$b=3$、$TOL=0.00001$ の場合には次のコードから $n=18$ になることが分かり、演習2の結果に一致します。
a = 0;
b = 3;
TOL = 0.00001;
log2((b-a)/TOL)-1 %MATLABでは、2を底とする対数を求める関数log2が用意されている
%他に自然対数を求める関数logと常用対数を求める関数log10があるが、
%底が任意の対数を求める関数は無いため、必要に応じて底の変換公式を利用する
ans = 17.195
二分法のアルゴリズムにおいて、計算量の観点から関数 $f$ の値の計算回数は最小限に抑えるべきです。$f(a)$ の値を格納する変数 $fa$ と $f(p)$ の値を格納する変数 $fp$ を新たに用いることで、演習2のコードを改良してください。
function value = f(x)
value = x.^3-2;
end
a = 0;
b = 3;
TOL = 0.00001;
fa = f(a); %関数fを使って、変数faの初期値を与える
n = 0;
while b-a >= 2*TOL
n = n+1;
p = (a+b)/2;
fp = f(p); %関数fを使って、変数fpの値を更新する
if fa*fp < 0 %関数fを使わずに、変数faとfpを使って条件判定をする
b = p;
else
a = p;
fa = fp; %関数fを使わずに、変数faの値を更新する(実はこの行は無くてもよい)
end
end
p = (a+b)/2;
printf("反復回数n=%d, 近似解p=%f\n",n,p)
反復回数n=18, 近似解p=1.259920
反復回数を $n$ とするとき、演習2のコードでは関数 $f$ を呼び出す回数は $2n$ 回ですが、演習4のコードでは $(n+1)$ 回に抑えることができています。 関数 $f$ が複雑であればあるほど、コード全体の計算量(実行時間)に差が生じます。