{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 二分法による方程式の近似解を求める\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 概要\n",
"\n",
"二分法によって、区間$[a,b]$における関数$f$の根の近似計算法を検討します。\n",
"\n",
"## 1. 二分法のアィディア\n",
"\n",
"関数$f=x^3-2$を例にして、区間$[0,3]$の中で$f(x)=0$の解の計算方法を説明します。\n",
"\n",
"準備として、関数fを定義します。 MATLABでは、以下の文法で関数を定義します。\n",
"\n",
"```\n",
"function 戻り値の変数 = 関数名(入力変数)\n",
" %関数の中では、関数値を計算して、その値を「戻り値の変数」にします。\n",
"end\n",
"```\n",
"\n",
"以下の関数の例では、xが入力変数、valueが戻り値の変数（出力変数）です。\n",
"\n",
"```\n",
"function value = f(x)\n",
" value = x.^3 -2\n",
"end\n",
"\n",
"```\n",
"\n",
" 補足: 「x.^3」の中の「.」演算子の役割\n",
"\n",
" xが配列であるとき、「x.^3」はxの各成分の3乗を計算して、計算結果を配列として返します。ここで、「.」という演算子によって、配列の成分毎の計算ができます。\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"%関数ｆの定義\n",
"function value = f(x)\n",
" value = x.^3 -2;\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"以下、$y=f(x)$のグラフを描いてみます。\n",
"\n",
"\n",
" 補足: 「hold on」の役割\n",
" \n",
" 二つのグラフを一つの軸に描くために、「hold on」という命令を使用します。以下の例では、y=0とy=f(x)を描いています。\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p = 1.5000\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot -f \"svg\"\n",
"%グラフ作成のフォマットを指定します。\n",
"\n",
"a=0;\n",
"b=3;\n",
"x=a:0.1:b;\n",
"y=f(x);\n",
"\n",
"hold on\n",
"grid on\n",
"plot(x,y*0,'b-') % x軸を描く\n",
"plot(x,y,'r-') % y=f(x)のグラフを描く\n",
"\n",
"p=(a+b)/2\n",
"plot(p,f(p),'r+') % pにおける関数の値を描く"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## グラフによる考察：\n",
"- 区間の両端では、$a=0,b=3$, $f(a)<0$, $f(b)>0$ が分かります。中間値の定理によって、$a$と$b$の間に$f(x)=0$の解が少なくでも一つが存在します。\n",
"- $p=(a+b)/2$とすると、$f(p)>0$が分かります。\n",
"\n",
"以下のコードでは、$a$, $b$を更新して、より小さい区間で近似解を検討します。\n",
"\n",
"- $f(a)$と$f(p)$の符号が異なるなので、$a$と$p$の間に$f(x)=0$の解が存在します。この場合、$b$を$p$にすると、新しい区間$[a,b]$で解の検討ができます。\n",
"- 2分法における重要な判定条件式：\n",
"\n",
"$$\n",
"f(a)\\cdot f(p)<0\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p = 1.1250\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot -f svg\n",
"a=0.75;\n",
"b=1.5;\n",
"x=a:0.001:b;\n",
"y=f(x);\n",
"\n",
"hold off\n",
"plot(x,y,'r-')\n",
"hold on\n",
"plot(x,y*0,'b-')\n",
"\n",
"p=(a+b)/2;\n",
"plot(p,f(p),'ro')\n",
"grid on\n",
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上記の計算を繰り返しすると、解の存在範囲を小さくすることができます。$a$と$b$が十分近いであれば、$f(x)=0$の近似解を$(a+b)/2$にすることができます。\n",
"\n",
"### 演習１\n",
"\n",
"$|b-a|<0.02$まで、上記のコードを手で繰り返して計算してください。\n",
"ただし、毎回の計算では、aを更新するか、bを更新するか、手動でコードを調整することが必要です。\n",
"途中の計算結果を以下のセルに書いてください。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. a=0, b=3, p=1.5, f(a)*f(p)<0\n",
"2. a=0, b=1.5, p=0.75, f(a)*f(p)>0\n",
"3. ...\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. プログラムで二分法のアルゴリズムを書く"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上記の操作をwhile文のループで行うと、効率な計算ができます。\n",
"\n",
"\n",
"終了条件:\n",
"\n",
"- [1)] 許容誤差$TOL$: $(b-a) N$ ($i$が計算の回数を数える)\n",
"\n",
"\n",
"\n",
"**計算アルゴリズム**\n",
"\n",
"\n",
"- 入力： 区間の両端$a$と$b$;許容誤差$TOL$; 最大計算回数$N$。\n",
"- 出力： 近似解$p$或いはエラー情報\n",
"\n",
"アルゴリズムの計算の流れ\n",
"\n",
"- [Step 1:] a, bを設定します。$i=1$, $FA=f(a)$ （初期化）\n",
"- [Step 2:] While ($i\\le N$ AND $(b-a) >= TOL$) do Step 3-6\n",
"- [Step 3:] Set $p=a+(b-a)/2$, $FP=f(p)$.\n",
"- [Step 4:] If $FP==0$ Then break the while. (Whileループを終了する)\n",
"- [Step 5:] If $FA\\cdot FP>0$, Then set $a=p$, $FA=FP$; Else set $b=p$.\n",
"- [Step 6:] Set $i=i+1$.\n",
"\n",
"- [Step 7:] If $(b-a)= TOL )\n",
" ....\n",
"end\n",
"\n",
"if (b-a)