Skip to content

N
numer-60-2
Commit f7f9d76a authored by Littichai Buddaken's avatar Littichai Buddaken
Browse files
Options

add week12(edited)&week13

parent a67cb5bd
Hide whitespace changes
Inline Side-by-side
Showing with 1417 additions and 0 deletions
week12/.ipynb_checkpoints/Numerical Integration-checkpoint.ipynb 0 → 100644
View file @ f7f9d76a
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Numerical Integration\n",
"\n",
"> คำนวณค่า $\\int_{a}^{b}f(x)dx$ เมื่อกำหนดฟังก์ชัน $f(x)$ มาให้\n",
"\n",
"คำนวณหาพื้นที่ใต้กราฟ $f(x)$ จากในช่วง $x \\in [a,b]$ โดยปกติแล้วจะใช้วิธีประมาณค่าโดยการแบ่งพื้นที่ออกเป็นช่วงๆ แต่ละช่วงจะเป็นสี่เหลี่ยมแล้วรวมพื้นที่สี่เหลี่ยมเหล่านั้นมาประมาณค่าพื้นที่ใต้กราฟหรือค่า $\\int_{a}^{b}f(x)dx$ นั่นเอง\n",
"\n",
"1. แบ่งเป็นสี่เหลี่ยมที่มีความกว้างเท่ากัน: Newton-Cotes\n",
"\n",
" 1.1. Trapezoidal rule\n",
" \n",
" 1.2. Simpson's rule\n",
" \n",
"2. แบ่งตามความเหมาะสม: Gaussian\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton-Cotes\n",
"$$\n",
"I = \\sum_{i=0}^{n}A_if(x_i)\n",
"$$\n",
"where\n",
"$$\n",
"Ai = \\int_{a}^{b}l_i(x)dx, i=0,1,2,...,n \\\\\n",
"x_i \\in [a,b]\n",
"$$\n",
"\n",
"### Trapezoidal rule\n",
"$$\n",
"I_i = [f(x_i) + f(x_{i+1})]\\frac{h}{2} \\\\\n",
"\\int_{a}^{b}f(x)dx \\approx I = \\sum_{i=0}^{n-1}I_i = [f(x_0) + 2f(x_1) + 2f(x_2) + ... + 2f(x_{n-1})+f(x_n)]\\frac{h}{2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"1. จงเขียนโปรแกรมเพื่อหา $\\int_{0}^{3}f(x)dx$ ด้วยวิธี trapezoidal rule โดยกำหนดให้ $f(x) = 3x^2 + 9x + 2$ โดยใช้ $n=4$\n",
"\n",
"```python\n",
"def f(x): return 3*x**2 + 9*x + 2\n",
"def trapezoid(f, a, b):\n",
" pass\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"74.34375"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def f(x): return 3*x**2 + 9*x + 2\n",
"def trapezoid(f, a, b):\n",
" n = 4\n",
" h = (b-a)/n\n",
" #sumI = 0.0\n",
" #for i in [0,1,2,3]:\n",
" # sumI += (1/2)*h*(f(a+i*h) + f(a+(i+1)*h))\n",
" return sum([ (1/2)*h*(f(a+i*h)+f(a+(i+1)*h)) for i in range(4) ]) \n",
"\n",
"trapezoid(f, 0, 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simpson's rule\n",
"* Simpson's 1/3 rule (First Simpson's Rule)\n",
"ประมาณค่าพื้นที่ใต้กราฟในแต่ละช่วงด้วย Larange Parabolic Interpolation โดยแต่ละช่วงจะเพิ่มจุดกึ่งกลางเข้าไป\n",
"\n",
"> จาก Larange Interpolation\n",
"> $$\n",
" P_n(x) = \\sum_{i=0}^{n}{y_i l_i (x)} \\\\\n",
" l_i(x) = \\frac{x-x_0}{x_i-x_0}\\cdot\\frac{x-x_1}{x_i-x_1}\\cdots\\frac{x-x_{i-1}}{x_i-x_{i-1}}\\cdot\\frac{x-x_{i+1}} {x_i-x_{i+1}}\\cdots\\frac{x-x_n}{x_i-x_n} \n",
"$$\n",
"\n",
"ถ้าใช้ Larange Interpolation จาก $[a,b]$ จะได้\n",
"$ x_0 = a, x_1 = \\frac{a+b}{2}, x_2 = b $\n",
"เมื่อทำการแทนสมการแล้วลดรูปจะได้\n",
"$$\n",
"l_0(x) = \\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} \\\\\n",
"l_1(x) = \\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} \\\\\n",
"l_2(x) = \\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \n",
"$$\n",
"พื้นที่ใต้กราฟในแต่ละช่วงคือ\n",
"$$\n",
"A_i = \\int_{a}^{b}l_i(x) = \\int_{-h}^{h}l_i(z)dz\n",
"$$\n",
"เมื่อกำหนดให้ $z_0 = -h, z_1 = 0, z_2 = h$ จะได้\n",
"$$\n",
"$$\n",
"นั่นคือ\n",
"$$\n",
"A_0 = \\int_{-h}^{h}\\frac{(z-0)(z-h)}{(-h)(-2h)}dz = \\frac{1}{2h^2}\\int_{-h}^{h}(z^2-hz)dz = \\frac{h}{3} \\\\\n",
"A_1 = \\frac{4h}{3} \\\\\n",
"A_2 = \\frac{h}{3} \n",
"$$\n",
"> คะแนนพิเศษ derive $A_1, A_2$\n",
"\n",
"$$\n",
"I = \\sum_{i=0}^{2}A_if(x_i) = [f(a) + 4f(\\frac{a+b}{2}) + f(b)]\\frac{h}{3}\n",
"$$\n",
"ดังนั้น first Simpson's rule สำหรับ $n$ ช่วงจึงเขียนได้เป็น\n",
"$$\n",
"I_i = [f(x_i) + 4f(\\frac{x_i+x_{i+1}}{2}) + f(x_{i+1})]\\frac{h}{3} \\\\\n",
"\\int_{a}^{b}f(x)dx \\approx I = [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... + 2f(x_{n-2}) + 4f(x_{n-1})+f(x_n)]\\frac{h}{3}\n",
"$$\n",
"\n",
"เงื่อนไขของการใช้สมการข้างต้นคือ $n$ ต้องเป็นเลขคู่ ถ้าเป็นเลขคี่แนะนำให้ Simpson's 3/8 rule [Second Simpson's Rule] กับ $x_0, x_1, x_2, x_3$ หรือ $x_{n-3},x_{n-2},x_{n-1},x_n$ ตามสมการต่อไปนี้\n",
"$$\n",
"I = [f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]\\frac{3h}{8}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"n,a,b = 4,0,3 \\\\\n",
"h = \\frac{b-a}{n} \\\\\n",
"x_i = x_0 + ih \\\\\n",
"\\int_{a}^{b}f(x)dx = [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3)+f(x_4)]\\frac{h}{3}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"1. จงเขียนโปรแกรมเพื่อหา $\\int_{0}^{3}f(x)dx$ ด้วยวิธี Simpson's rule โดยกำหนดให้ $f(x) = 3x^2 + 9x + 2$ และ $n=4$\n",
"\n",
"ตอบ:\n",
"\n",
"```python\n",
"def f(x): return 3*x**2 + 9*x + 2\n",
"def simpson(f, a, b):\n",
" a,b = 0,3\n",
" n = 4\n",
" x0 = a\n",
" h = (b-a)/n\n",
" return (f(x0+0*h) + 4*f(x0+1*h) + 2*f(x0+2*h) + 4*f(x0+3*h) + f(x0+4*h))*h/3\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"73.5"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def f(x): return 3*x**2 + 9*x + 2\n",
"def simpson(f, a, b):\n",
" n = 4\n",
" x0 = a\n",
" h = (b-a)/n\n",
" return (f(x0+0*h) + 4*f(x0+1*h) + 2*f(x0+2*h) + 4*f(x0+3*h) + f(x0+4*h))*h/3\n",
"simpson(f, 0, 3)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def simpson(f, a, b, n):\n",
" h = (b-a)/n\n",
" x0 = a # x_i = x_0 + i*h\n",
" #f(x0) + 4*f(x1) + 2*f(x2) + .... + 4*f(x_n-1) + f(x_n)\n",
" s = f(x0)\n",
" \n",
" # odd - x1, x3, x5, ... \n",
" #s += 4*f(x1) + 4*f(x3) + 4*f(x5) + ...\n",
" for i in range(1,n,2):\n",
" s += 4*f(x0+i*h)\n",
" \n",
" # even - x2, x4, x6, ...\n",
" #s += 2*f(x2) + 2*f(x4) + 2*f(x6) + ...\n",
" for i in range(2,n-1,2):\n",
" s += 2*f(x0+i*h)\n",
" \n",
" s += f(x0+n*h)\n",
" return s*h/3"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(4, 0.6321341753205324),\n",
" (6, 0.6321232596014171),\n",
" (8, 0.6321214146047422),\n",
" (10, 0.6321209095890152),\n",
" (12, 0.6321207280453661),\n",
" (14, 0.6321206501876014),\n",
" (16, 0.6321206123891728),\n",
" (18, 0.6321205922694484),\n",
" (20, 0.6321205807706575)]"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from math import e\n",
"def f(x): return e**(-x)\n",
"a,b = 0,1\n",
"simpson(f, a, b, 4)\n",
"[ (n,simpson(f,a,b,n)) for n in range(4, 22, 2) ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercise\n",
"\n",
"1. จงเขียนฟังก์ชัน `def simpson(f, a, b, n)` เพื่อคำนวณหาพื้นที่ใต้กราฟของฟังก์ชัน $f$ จากค่า $x=a$ ถึง $x=b$ โดยแบ่งเป็น $n$ ช่วง ดัวยวิธี Simpson's rule\n",
"\n",
"2. จงเขียนฟังก์ชัน `def trapezoid(f, a, b, n)` เพื่อคำนวณหาพื้นที่ใต้กราฟของฟังก์ชัน $f$ จากค่า $x=a$ ถึง $x=b$ โดยแบ่งเป็น $n$ ช่วง ด้วยวิธี trapezoidal rule\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import sympy\n",
"sympy.init_printing()"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def f(x): return 3*x**2 + 9*x + 2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2, 56)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(0), f(3)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from math import e\n",
"z = e**-x"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAALYAAAARBAMAAACCx0thAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADGUlEQVQ4EbWUz2scZRjHP5Pdmc3uZuMa8KCC\nrlETQaSVGA+ehoIeRMicKmohy6LF1pgsVRJQ6e5FBD3EUgv+OHRFUMFDFqT1IiQN1h4UunooeLHr\npaUQtDVGjUHH7/u8E91/IM9h5nm/73c+7zPvO8/A3kRh/Yln94YM56PupwTjU7Fb4PHG0lLTVgq/\nPdaGsbsfgkMX7jMJTDSvn4YD8GVyoEc4uVjV5WpM+Xz+gcxOXO7DQQp/OuGXNE3rVE7BELxIWKfV\npkbUgx+STDSvTeuBj+B0+gc81Q9OUO6FH1JM0yaVKcU03FSF5+A9x74LIg5ObMGPcDtRQvHtSp+g\nGzYWxDbRvJbByEn4Yk6AI9DgGrxKefGcI/mYIeZ9WIg1bMLLekLsl+Bpih2i3wr7CTtwUWwTzWsZ\n3PE51FHsQIuH4XISuaGP4PorIq4mng3Da569PM1xyltiD28mhb5nm2hey1RLxg43ZYh/h+XeADuc\ne+OCKoIVu5JT6uqO0qt9pZS2WPjnuBJX9664kvhstCr2/MY98IHqfv1Xsdein9dj2QdjWAfiwmEc\nm9aOLbbcpJCe0dixM9F5bfoQYk9zZ1sbwupPkmeauSSwz0Lu3SjVLQuv6+bYledX3nXKx3Dpze22\nZ2eivD5rOjbkOpTXKq3vVPdMU+NHnDgQ4z4f6eru2FcY2ValQzVG61w85dmZKK9lldizh27Axlet\nR21P9Py+6gBYuJoflvbr7tgnVEMP1qFYJf+3Z3vReS17ErGjbtYbC4nO8nLvBe1X28Oy6/cEtths\nXYLYbm+iPpUaN89KmTd2Jsp7i01fOnx4e67YZUi1wFmm1HyJTnVf4sa7ka8xYuyZmqSs7lzVtfQ5\n1c2D/9Wdq5pXdWsaTroGLnW5Vs3fcL1zv7qYzzTxf9y2dOwos129U02iY7+VsEHwydJEZ/gMo2ue\nbaJ5LZN1k6AmMN8kz/Qp9cJ3GGNUCw/Eapr+hftUbu3p/I7szJNv6F9V1u+lw2OT+ldNnD7b96J5\nbVqNms4ztngvFMblCV/7OiaYbCQD5D1N/wX8wRK5TpgTOQAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$2.71828182845905^{- x}$$"
],
"text/plain": [
" -x\n",
"2.71828182845905 "
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAARBAMAAAARYzyGAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMlTvq5l2Zoki\nu0Rn3bgMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADUElEQVRIDcWVS0wTURSG/yktQ9spNLpA4qIl\nMTFoDES3Kk1ciHHRrsToAl9xYVDqE2iMVKKBBMNjJQQTutIgCxs1ASHGURNC3Ahx4UaSRpcueAho\nQBjPOXc6Q4C1nMXMud898//3zJ3eAtsQwd4Tt7bBFujRo0c3GWsdQzEHGu3TMjjWNDaWFar1jeSg\ncFl/p13oQsmI1gAV4ZoEdo3dBUTS0+N9aJcj5snkU7r7zsvAk9AOO7QUwTUevLAsKykVRcBtKJyG\nngBOhQEXSkYPHAEarBVoVzBuQiQLLSsL3zTFFFBiUo0dpW1zkj0FXuYZrgOXefAI0CEVZ4APCvsy\nMKJaUwsZu1AywH8V2N9sQj8HPQ2R9Aw/cWQRR8wd+JXxV6COtFRcBFpilGaB1yRHFW+BnRAcnIRW\nBQxQsQslAz4dAJL0XEkUwVmIpK4U+WrMvmLFfNjGv4GZRJ5VhpUxUNytjGemMArBxavhYEYZu1Ay\nWqhtHJmEfwEiuc5Ya/7ZH85bKFn6EP6QMZk4Ua1KQgR4abo1lOE5wi3Lo5Rwxy5UWcAk49TgAxRS\nx2tKUn/eG+MHN4fq2FgB4ll3tpjGHOwhFe/meSWMg1Y5ZWwMF0pWBjKewuccvWfPspIMhQ35Tql6\nQ9jG1LEYa2q6ICl3bZZuXOG7V32NUsa1P5ZyytiFKsuyMRCqQgXOLhp5yS8MndhxiGPPplcdiUpJ\nuyr084iNn8G/RD0SDiQxQL9A7tiFkvliyrjoLwJ9Jxec3esyldiGq+qYv4S6BE2pVv1pVVUwSXeu\nqKcXkgDjQhPeRWXsQsl2g4x13l56SJ9TkjdoQ3I03hy28TQdF9SEHadhmJxGknShCn7legaMI0RS\nYuzC9zJd29i41EzfVREtFJ4kRPIC0OXq0oQTtjH92lsd5k3Db/IonqaL3XHI9DIuNIFvTsc2pI5D\nxHGVT7OCaKAecZMPkFakgYM0sUWwLG1sQUK75Mx+HBu5yRB1aWJc8T2MQQguLkegWxm7UDIqXYVB\nJ5bpL9f20QJYcgK0iq2i6M58ijdWe/M45sxXWtYv2e3xBJ3IXOFtoj8JwTje0Qm0NezNrIMyTUeZ\nlcLE8H1gpD1HZwNLGh1NW79px+w/Jf8Ai2Q6WOYdK4QAAAAASUVORK5CYII=\n",
"text/latex": [
"$$- 1.0 \\cdot 2.71828182845905^{- x}$$"
],
"text/plain": [
" -x\n",
"-1.0⋅2.71828182845905 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sympy.diff(z)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
week12/Numerical Integration.ipynb 0 → 100644
View file @ f7f9d76a
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Numerical Integration\n",
"\n",
"> คำนวณค่า $\\int_{a}^{b}f(x)dx$ เมื่อกำหนดฟังก์ชัน $f(x)$ มาให้\n",
"\n",
"คำนวณหาพื้นที่ใต้กราฟ $f(x)$ จากในช่วง $x \\in [a,b]$ โดยปกติแล้วจะใช้วิธีประมาณค่าโดยการแบ่งพื้นที่ออกเป็นช่วงๆ แต่ละช่วงจะเป็นสี่เหลี่ยมแล้วรวมพื้นที่สี่เหลี่ยมเหล่านั้นมาประมาณค่าพื้นที่ใต้กราฟหรือค่า $\\int_{a}^{b}f(x)dx$ นั่นเอง\n",
"\n",
"1. แบ่งเป็นสี่เหลี่ยมที่มีความกว้างเท่ากัน: Newton-Cotes\n",
"\n",
" 1.1. Trapezoidal rule\n",
" \n",
" 1.2. Simpson's rule\n",
" \n",
"2. แบ่งตามความเหมาะสม: Gaussian\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton-Cotes\n",
"$$\n",
"I = \\sum_{i=0}^{n}A_if(x_i)\n",
"$$\n",
"where\n",
"$$\n",
"Ai = \\int_{a}^{b}l_i(x)dx, i=0,1,2,...,n \\\\\n",
"x_i \\in [a,b]\n",
"$$\n",
"\n",
"### Trapezoidal rule\n",
"$$\n",
"I_i = [f(x_i) + f(x_{i+1})]\\frac{h}{2} \\\\\n",
"\\int_{a}^{b}f(x)dx \\approx I = \\sum_{i=0}^{n-1}I_i = [f(x_0) + 2f(x_1) + 2f(x_2) + ... + 2f(x_{n-1})+f(x_n)]\\frac{h}{2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"1. จงเขียนโปรแกรมเพื่อหา $\\int_{0}^{3}f(x)dx$ ด้วยวิธี trapezoidal rule โดยกำหนดให้ $f(x) = 3x^2 + 9x + 2$ โดยใช้ $n=4$\n",
"\n",
"```python\n",
"def f(x): return 3*x**2 + 9*x + 2\n",
"def trapezoid(f, a, b):\n",
" pass\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# โจทย์ f(x)= e^x*sin(X) ,n = 10\n",
"def f(x) :\n",
" from math import e ,sin\n",
" return e**x*sin(x)\n",
"#n คือ จำนวนการแบ่งช่วง\n",
"#h คือ ค่าของช่วงที่แบ่ง from (b-a)/n\n",
"def area1(f,a,h,n):#Findin' area from f , a , h , n \n",
" return sum([(1/2)*h*(f(a+i*h)+f(a+(i+1)*h)) for i in range(n)])\n",
"def area2(f,a,b,n):#Findin' area from f , a , b , n without h\n",
" return sum([(1/2)*((b-a)/n)*(f(a+i*((b-a)/n))+f(a+(i+1)*((b-a)/n))) for i in range(n)])\n",
"area1(f,0.1,0.05,10),area2(f,0.1,0.6,10)\n",
"#Chernaling are you kiddin' me\n",
"#Loadto be the s'one\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.2574269780305713 0.2574269780305713\n"
]
}
],
"source": [
"import sympy\n",
"sympy.init_printing()\n",
"a, b =0.1,0.6\n",
"n, h =10,(b-a)/10\n",
"def y(x):return 3*x**2+9*x+2\n",
"print(area1(f,a,h,n),area2(f,a,b,n))"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAAVCAYAAADo8PLxAAAABHNCSVQICAgIfAhkiAAACv1JREFU\neJztnX2QVlUdxz8ICGS2FL2QSa3hG4yNaIlRps9m4WBa9jLmH02uZkWk5RDaQGNub0SWhlkN2FRE\nZTkDg4wjIYiKLoZByWQIWcCDBquJDC8VoUvbH99z57nP3Xvvc+65d5f77J7PzM599rzd83vuOd/n\nnt85517weDwej8fj8Xg8Hk+/8Qvgn8BxOct5O9ADXJO7Rh6Px1MsrUifFhVYZtNp58eAO4BHgQPm\npL9yLOs7wBrgWeAQsBd4ArgZGJOQp2rOGff3XEKeIcCngceBfwH/BjYC04FjYtIXYeMnQvVKuygn\nAj8DdgOHkX3zgVcnpC+jLWNM+DLg7+ha7gc6gU8l1MslT1bbIXt7calXmTgH+B8wMybOpS0sA7qA\nVxZYx6xk7SNJVMmuHX3NYNBTyG+n11Ovp1FaKfaGrCm1c5OpzEFgi2XFkngJWI860Dxk8AZT5i5g\nXEyeKrAP6Ij5m5Vwnl+bMp8HfgLcDjxlwhbHpM9r4zhTx4OkC8h4U6ce4B70HTxo/t9KvIiW0Zbp\nJm63qd+30TXdZ8KXoM6fN09W2yF7e3GpV5lYheo6KibOpS1MNunmFFjHLLj0kSSqZNeOvmYw6Cnk\ns9PrqdfTOIYDpwNvdMwfpSm1sw04BX2JFcuKJTEyIfxbptwfx8RVzZ8tHzZlbQdeGwo/FrjXxH0k\nkiePjUOAB4BtwHdJF5D7Tfx1kfDbTPiCJrHlvcCl9B7tjAWeMfk+mjOPi+2Qvb242FIWTkUjvDsT\n4l3bwhZgJ26j2XZznopDXsjeR9Kokq0t2NKOu42DQU/B3U6vpzW8nvYdZdTOzFTIJyBJnGnKXR0T\nVyVbg1hsyvp8TNwkE/dgSv4K2Wz8Irqw56NRQ1KnG2/idtD7Yh1PzYUcnssuqy1pzDH57siZx9X2\nKsX9CKfZMhm4G3kiDiNX9Srg8pi0lwOPINf9IeBJYDYwIuG8H0RTUV2m7N3AWmBGJN08U78LLWyp\nYN8WbjZpL7JIG6Ud95sVlz6SRpXy3ZCFqTDw9RSy2en1tJ7Boqc2mtdq8i+K5A2HtwK/BfYA/0VT\nspfEnK/ftLPsa17iuNQc/5wQPwLNw89BjbwNGJqQdqw5bo+JC8Leg0YGeZmALuzt6Ac3jTZzXIU6\naZiDwDrgFcA7Q+FltSWNl82xO2eePLZnaS9Z6wVah/EYcJk53grcB7ye3jdNc9GN2wTgLuCHaNQ1\nF43wo/X/DLAcmIhGrrcCK5Bb/apI2vcBR9C0VZGsM8f3F1xuI1z6SCOKagvNhNdTr6fRPGXW0yya\nl8ZbgD+gm7JfIt09w5TdFknbb9o5rOAT9AWz0MK3FuAdwHlIPOYlpB+LvuAwO9DFWhsJ32OOJ8WU\n81ZzHGY+b81U63qGmTo9g92c8Wnm+HRC/N+AqciVusaEldWWtHI+aT6vzJknj+1Z2kvWek1EU0EH\nkIBtjuQ7MfR5CvKEPYs8asEi2NloAeglqC/MDeX5LFobdCba/RMmPNVwHBrZbkGegCLZYI7nF1xu\nI1z6SCOKaAtlx+tpb7yeNoeegr3mNaKCPJFfC4XdZc53A/CQCetX7WwGD9ks5Nq7HonHStR5XohJ\n+3PkVhyLvsi3AQvRXfDv0EUMc585zgReEwofTv2FyrpjK8pXgbPQ9MUhi/Qt5rg/IT4IHx0KK6st\nScxDI5IVyPuTJ4+r7VnbS9Z6fQ6JyzfofTMG8I/Q56vN8ZvU70jqBr6ERvZxUxjd1EaTYfaEPr8J\njVK7Ei1wZz9y97+5D8pOw6WPpFFUWyg7Xk974/W0OfQ0wEbzGrETaW2Y+9GN8eRQ2FHVziqa07T9\nS5snrVikycIb0GLDv6I547Mz5P2eqcuySPhQJEg96EdwIXIbb0bbwneauHMTyq3Q2MZzUQO6JRLe\nQfI6gTtT4qC2EHd2E9gSxxdM+i3Ud3jXPHltj5LUXrLWK9jFdrpFOX80aU9OiA9saAmFzTRhXcD3\n0bTo62LyTjHp7raoB2Tvu7toPE1SJZu2LGpQnksfcSFLW6hSrI1hKgx8PYXGdno97c1g0VOw17xW\n0teQ3ZNw/k40PRnQr9oZnbLchu7YbNmdIW1enkcX9E/I9bwY3UXbsAB5GaLTKkfQGoqZaN77SmT/\nw2h3xxKTLuoatWWYqefTwE0Z8gUjtpaE+CB8XyisrLZEuZbaFuoLUQfPm6do25PaS9Z6BSPuXRbn\nDK5p0kisC42kRlNrH7ehUeEMJGTXIzFYi9zuG026YOSdtNMuL6NoPLqfT29v1STgQ+iBi9VI3KYG\n5bn0ERds2wIUb2Nf4vW0htfTesqqp2CveY1I0oVu6mcOy6CdhVChb3YFgR5o2IP9nHGLSZ/lZnMk\n2sER58oPqJBu42jsR8vzQ/muMWELE8oNtnDb7PqAo2tLmKDzPIkWtdvgkieMje1RbNqLTb1cPGTj\nE+KDUWnSFNxo4GL0zKAjwIvURo4nmLydFvWAbH33GDSdus2y7DDt5jwVh7xF95EkXLQjTDvuNoap\nMPD1FNLt9Hpaz2DT0yhpmtdKuocsGh7wsIkP6FftbIZF/XGcYI5HUlPVCHbPxO0aSeIKtIvkNxny\nRDkM/DQh7my0dqATTRv8PhT3kDlOpXbRAo4H3g38B/tdH0fTloAvo7UBm9CuEpv5fpc8UVxsb9Re\nbOu1Hi2cnkbjhb9PoO+xQu+bm5PRBoAdJI/s9qE1FytQm7kajUiXIu/aC9QWNxfJaWgnaH97e4ru\nI0m4aEez4fXU66ktR1NPo6RpXlEMGO2sYHenOB55EIaHwk4l3r18DLX5/nWRuAnEP3OoFe2i6SF+\nF8urYsImoYuwl5pYxVHBfdTaQfo6AZeHXpbVlptM/Ebs1zhkzZPVdtf2kqVeE9Hi073mc5TwLst3\nmXJ3UL8mYiha79ADfCWSv434J1kHD2+cFgpbYsKS1qiFqWDfFq4yaa+1SBulnXzeI5c+Eqc3rm3B\nhnb610PWzHoK7jrUgdfTNAaCnoK95rVSjIcM+lE7i/aQXWb+oPYskynUjN9D71corEHPBDmJ2hqL\ni9FrFDrRD9SLaBHqBWir7XPo+U5hPo7mqR9B0zsHkTh9ALlZV6DFhVFWozncv5g8E0yeQ2gePbpO\nzsVGF2ag51b9ALnSt6BFlG1o3UH0xxnKacuVwNfR6PtRNO8fpUp9B3HJk9V2l/aStV5Poeu4AHnA\nliNxGoPejXaA2jNvHkOLe280NixB26ynobU9negJ3mGWoYdarjfnHYIer3EOmgJ9IJR2KVr/cRF6\nb1wU17YwFX0fy2Pi+hqXPhKnN67a0dcMBj11tTMrXk+bX08hm+YVRdNqZwfpc+HVmDxVE9caCjsD\nPRRzEzKmGy3M3GDOEXcnfQFyo25FrsyX0d38avRMk6R3Yt2ALuQ+5EbeDvyIeu9FXhvTyknbSTMO\nbSXuQs9e2Un6y3DLaEujc/SgUUnePFltd2kvLvUCdc6laBHsS0jMVqKX0ka5Av1wHkRrLjajH4u4\nRaXTkUBtR1MuwQujb0RTMWGORQu5H48px8a2akyeFiTQSTuWGtFOfu9R1j5SpbfeuGqHDe2429jB\nwNdTVzuTyvB6OrD11FbzWinOQ1ZG7fR4PE3ObCQQZxVU3nWmvPMKKs/j8XjKiNdOj8dTKCORV+De\nAsoahTx9Sxol9Hg8nianX7RzoL+nzePx1OhGz50agRbRxj3t2pZT0FTGLeR/1pfH4/GUGa+dHo/H\n4/F4PB6Px+PxeDweT5/zf8ACgG8g7vYxAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$- 1.35914091422952 \\cos{\\left (1 \\right )} + 0.5 + 1.35914091422952 \\sin{\\left (1 \\right )}$$"
],
"text/plain": [
"-1.35914091422952⋅cos(1) + 0.5 + 1.35914091422952⋅sin(1)"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x,y= sympy.Symbol('x'),sympy.Symbol('y')\n",
"y = 3*x**2+9*x+2\n",
"sympy.integrate(y,(x,0.1,0.6))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.79390625"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def f(x): return 3*x**2 + 9*x + 2\n",
"def trapezoid(f, a, b):\n",
" n = 4\n",
" h = (b-a)/n\n",
" #sumI = 0.0\n",
" #for i in [0,1,2,3]:\n",
" # sumI += (1/2)*h*(f(a+i*h) + f(a+(i+1)*h))\n",
" return sum([ (1/2)*h*(f(a+i*h)+f(a+(i+1)*h)) for i in range(4) ]) \n",
"\n",
"trapezoid(f, 0, 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simpson's rule\n",
"* Simpson's 1/3 rule (First Simpson's Rule)\n",
"ประมาณค่าพื้นที่ใต้กราฟในแต่ละช่วงด้วย Larange Parabolic Interpolation โดยแต่ละช่วงจะเพิ่มจุดกึ่งกลางเข้าไป\n",
"\n",
"> จาก Larange Interpolation\n",
"> $$\n",
" P_n(x) = \\sum_{i=0}^{n}{y_i l_i (x)} \\\\\n",
" l_i(x) = \\frac{x-x_0}{x_i-x_0}\\cdot\\frac{x-x_1}{x_i-x_1}\\cdots\\frac{x-x_{i-1}}{x_i-x_{i-1}}\\cdot\\frac{x-x_{i+1}} {x_i-x_{i+1}}\\cdots\\frac{x-x_n}{x_i-x_n} \n",
"$$\n",
"\n",
"ถ้าใช้ Larange Interpolation จาก $[a,b]$ จะได้\n",
"$ x_0 = a, x_1 = \\frac{a+b}{2}, x_2 = b $\n",
"เมื่อทำการแทนสมการแล้วลดรูปจะได้\n",
"$$\n",
"l_0(x) = \\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} \\\\\n",
"l_1(x) = \\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} \\\\\n",
"l_2(x) = \\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \n",
"$$\n",
"พื้นที่ใต้กราฟในแต่ละช่วงคือ\n",
"$$\n",
"A_i = \\int_{a}^{b}l_i(x) = \\int_{-h}^{h}l_i(z)dz\n",
"$$\n",
"เมื่อกำหนดให้ $z_0 = -h, z_1 = 0, z_2 = h$ จะได้\n",
"$$\n",
"$$\n",
"นั่นคือ\n",
"$$\n",
"A_0 = \\int_{-h}^{h}\\frac{(z-0)(z-h)}{(-h)(-2h)}dz = \\frac{1}{2h^2}\\int_{-h}^{h}(z^2-hz)dz = \\frac{h}{3} \\\\\n",
"A_1 = \\frac{4h}{3} \\\\\n",
"A_2 = \\frac{h}{3} \n",
"$$\n",
"> คะแนนพิเศษ derive $A_1, A_2$\n",
"\n",
"$$\n",
"I = \\sum_{i=0}^{2}A_if(x_i) = [f(a) + 4f(\\frac{a+b}{2}) + f(b)]\\frac{h}{3}\n",
"$$\n",
"ดังนั้น first Simpson's rule สำหรับ $n$ ช่วงจึงเขียนได้เป็น\n",
"$$\n",
"I_i = [f(x_i) + 4f(\\frac{x_i+x_{i+1}}{2}) + f(x_{i+1})]\\frac{h}{3} \\\\\n",
"\\int_{a}^{b}f(x)dx \\approx I = [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... + 2f(x_{n-2}) + 4f(x_{n-1})+f(x_n)]\\frac{h}{3}\n",
"$$\n",
"\n",
"เงื่อนไขของการใช้สมการข้างต้นคือ $n$ ต้องเป็นเลขคู่ ถ้าเป็นเลขคี่แนะนำให้ Simpson's 3/8 rule [Second Simpson's Rule] กับ $x_0, x_1, x_2, x_3$ หรือ $x_{n-3},x_{n-2},x_{n-1},x_n$ ตามสมการต่อไปนี้\n",
"$$\n",
"I = [f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]\\frac{3h}{8}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"n,a,b = 4,0,3 \\\\\n",
"h = \\frac{b-a}{n} \\\\\n",
"x_i = x_0 + ih \\\\\n",
"\\int_{a}^{b}f(x)dx = [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3)+f(x_4)]\\frac{h}{3}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"1. จงเขียนโปรแกรมเพื่อหา $\\int_{0}^{3}f(x)dx$ ด้วยวิธี Simpson's rule โดยกำหนดให้ $f(x) = 3x^2 + 9x + 2$ และ $n=4$\n",
"\n",
"ตอบ:\n",
"\n",
"```python\n",
"def f(x): return 3*x**2 + 9*x + 2\n",
"def simpson(f, a, b):\n",
" a,b = 0,3\n",
" n = 4\n",
" x0 = a\n",
" h = (b-a)/n\n",
" return (f(x0+0*h) + 4*f(x0+1*h) + 2*f(x0+2*h) + 4*f(x0+3*h) + f(x0+4*h))*h/3\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"73.5"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def f(x): return 3*x**2 + 9*x + 2\n",
"def simpson(f, a, b):\n",
" n = 4\n",
" x0 = a\n",
" h = (b-a)/n\n",
" return (f(x0+0*h) + 4*f(x0+1*h) + 2*f(x0+2*h) + 4*f(x0+3*h) + f(x0+4*h))*h/3\n",
"simpson(f, 0, 3)"
]
},
{
"cell_type": "code",
"execution_count": 97,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def simpson(f, a, b, n):\n",
" h = (b-a)/n\n",
" x0 = a # x_i = x_0 + i*h\n",
" #f(x0) + 4*f(x1) + 2*f(x2) + .... + 4*f(x_n-1) + f(x_n)\n",
" s = f(x0)\n",
" \n",
" # odd - x1, x3, x5, ... \n",
" #s += 4*f(x1) + 4*f(x3) + 4*f(x5) + ...\n",
" for i in range(1,n,2):\n",
" s += 4*f(x0+i*h)\n",
" \n",
" # even - x2, x4, x6, ...\n",
" #s += 2*f(x2) + 2*f(x4) + 2*f(x6) + ...\n",
" for i in range(2,n-1,2):\n",
" s += 2*f(x0+i*h)\n",
" \n",
" s += f(x0+n*h)\n",
" return s*h/3"
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 4 1.3616491974399736e-05 0.003288870199135485\n",
" 6 2.7007728591321722e-06 0.001462565055180165\n",
" 8 8.557761842498124e-07 0.0008228593819220587\n",
" 10 3.5076045723503313e-07 0.0005266793587330731\n",
" 12 1.6921680812131967e-07 0.0003657681764012155\n",
" 14 9.135904344148571e-08 0.0002687358901396486\n",
" 16 5.356061483219321e-08 0.00020575501594155554\n",
" 18 3.3440890390146194e-08 0.00016257408557796005\n",
" 20 2.1942099470706466e-08 0.0001316862962578158\n"
]
}
],
"source": [
"from math import e\n",
"from sympy import *\n",
"def f(x): return e**(-x)\n",
"a,b = 0,1\n",
"r = 0.632120558828558\n",
"simpson(f, a, b, 4)\n",
"for v in [ (n,simpson(f,a,b,n),area2(f,a,b,n)) for n in range(4, 22, 2) ] :\n",
" print('{:5} {:20} {:20}'.format(v[0],abs(v[1]-r),abs(r-v[2])))"
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAAASCAYAAADBs+vIAAAABHNCSVQICAgIfAhkiAAABgBJREFU\naIHt2l2sXUUVB/AfTbFYNXwJJQpaQIgIhoimohi5V7QEK6agGGIMPqiJJgoYFaKJoT4YUiSIwQga\noz5oSDQIJCix2vAlSpBIRRARU07V1KYU+Shf1tbjw5qTnu6z9zl75ux7H8z5JyeTu2f+a639nzuz\nZ9YMM8www4LhSHwPW/Fv9HA1Di60dwZuxLZkbyt+gfdW2q3HRvwdL+BfuB+X4dAau4fi48n2XxPn\nafwaH8OShng+iGtwF55BHz9saFvqg3wdeymWut+2DjlrsAH/SO+zGT/B28a8Sw5nGs1y/PQsjl6t\n49qvQjoWv8HhuBl/xirM4xGchifGOK3iCnwhBXErduAwvBm/wiVDbXfh9/gTtuNlOBVvEf+Qp4qB\nNsAncS3+idvwN6zAuTgQN+A8IdQwNuFkPJviej1+hI/UxF/qo0THHg4SA66KZ3FlzfNcznqh+RO4\nSfTH6/B+LMUFRieWXE6pZrl+ct+9lFOiGeKL0cdnKs+vSs+vqyM14BOJ8wO8pKZ+/8rfBzTY+Wqy\n863K83fhbKMz3RGiA/v4QI29eRwnJpI5479MpT5KdOylXw5yOEdgj5iBD6/Uzae4NnfAKdGsxE/P\nwupVGhdiNu3jMaNCvEKM3OfEF2MSlomvyxb1AykHJ6e4fpnB+VLiXDOh3Zzxg6nER6mOPQv7z/HW\nFNfNDfXPYGcHnHFo0qzET8/CD6asuJYOVcyncgP+WyHtxN1YLZZbGycE8R6xnLs62VqDk/Ai7sVv\nJ/CHcXYqH8jg/CeVuzM4uWjyMY2Oy8Ry8zViwD2AO8Xs2IS2nEfFUnoVXimWKwO8Uwz0mzrgjEOT\nZqV+FlKvaeLyNTEKP9cQxDdT/afGBDrAV1Lby/FHo5u9O8Rgq8PnsQ5fF0mCPv4wpn0VS4d8njmh\n7ZyyL9M4H6U69tRvjDfj9AZbuZyLxQDfju+I/vmxmOQ2GF3KlHLqMKlfcv30LLxeJXGRGvZFJqYO\ng73LFxvqh3FtartbjPx34OV4o737idsbuNvs+6K3ig1sW1yZeD9r0XZO2WAa56NUx8vEfmMFlosv\n+XWiI58Xy90qSjhrRaZ0WONH8eGGeEs5VbTplxw/i6VXblzodjB9O7V9ESsrdctFVq5vfDp2Bc4R\n2a+tOKWF3wuT3YdxSIv2c/IH0yQfXerI3n/CGzNibOJcIia4q3CM6ItT7J3grqixVcKpok2/dOGH\nbvUqjqvLZd761LZpb/TdVH9RC1uvFec0D05o9+lk8yGRhWmDOXmDqY2PLnUk0rB9eUcSdZy59Oyn\nNe2Xi2OCPeIfZhpOFW0068LPAF3plR3XcLbpkVQe3+DwuFT+pUVwA1tPNdQ/mcqXtrC1RZw9nSg2\ngXW4WGSIHhQJgHEHcKVo66NLHeHxVLbJoo7jvC+Vt9W0f14khpbgTVNyhtFWs2n9DKMrvbLjGh5M\nA8Jq9Snd05KBe1oEt1GM6DfU2CLWqkT6uA1elcq6jMulIlmxSXTY9pY2c5Djo0sdiawfDecZGZxl\nqWxK5Aye75qSM0COZtP4qaIrvaaOK/ew8Vhxg6B6AEvk5vv4bOX5arHhe1KciBOz+IFGscTePcbd\nNfVfTnX3abdHqmLO5GVeiY9cHU9QP5OuFBvdvjijmYbzofRsG15d4Zwl+uQF+17dKuGQr1mun8XQ\nKzuuSdeJHhYHV/NiWfJ2+64re2JPc7TRw7Ajk62jxJfq/tRubQrwfHG1hFgOXC7ubz2WfKwQ6cpj\n0sucIZZ7A3xU3K7YI5YSTxtFL7UZxtr0I9bwZ4oZ6a70bIdIz0/jI1fHdWKPdadY1u5MNtaImyE/\nF8mYXVNwlohB/u7UdnBf8gSxnNlP9MM3hnyUcEo0y/WzGHqVvv8+OArfF3erdiXHTRc0e2JgrGyw\ndZgQdEuytSMFtKrS7iSxMd+U2uwWnfC7JELd7LZO/ZnB8O/2Al6vAx/k6Xg6rhd3+J4Sh5uPi1sf\nFxid9Eo5+4vOv0ec3u8Wy69bxIqhDrmcdco0y/GzWHqVvP8MM8wwwwwzzDDDDP9f+B/5DlG+xc4w\nEgAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$0.632120558828558$$"
],
"text/plain": [
"0.632120558828558"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sympy import *\n",
"f = e**(-x)\n",
"integrate(f,(x,0,1))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercise\n",
"\n",
"1. จงเขียนฟังก์ชัน `def simpson(f, a, b, n)` เพื่อคำนวณหาพื้นที่ใต้กราฟของฟังก์ชัน $f$ จากค่า $x=a$ ถึง $x=b$ โดยแบ่งเป็น $n$ ช่วง ดัวยวิธี Simpson's rule\n",
"\n",
"2. จงเขียนฟังก์ชัน `def trapezoid(f, a, b, n)` เพื่อคำนวณหาพื้นที่ใต้กราฟของฟังก์ชัน $f$ จากค่า $x=a$ ถึง $x=b$ โดยแบ่งเป็น $n$ ช่วง ด้วยวิธี trapezoidal rule\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import sympy\n",
"sympy.init_printing()"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def f(x): return 3*x**2 + 9*x + 2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2, 56)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(0), f(3)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from math import e\n",
"z = e**-x"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAALYAAAARBAMAAACCx0thAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADGUlEQVQ4EbWUz2scZRjHP5Pdmc3uZuMa8KCC\nrlETQaSVGA+ehoIeRMicKmohy6LF1pgsVRJQ6e5FBD3EUgv+OHRFUMFDFqT1IiQN1h4UunooeLHr\npaUQtDVGjUHH7/u8E91/IM9h5nm/73c+7zPvO8/A3kRh/Yln94YM56PupwTjU7Fb4PHG0lLTVgq/\nPdaGsbsfgkMX7jMJTDSvn4YD8GVyoEc4uVjV5WpM+Xz+gcxOXO7DQQp/OuGXNE3rVE7BELxIWKfV\npkbUgx+STDSvTeuBj+B0+gc81Q9OUO6FH1JM0yaVKcU03FSF5+A9x74LIg5ObMGPcDtRQvHtSp+g\nGzYWxDbRvJbByEn4Yk6AI9DgGrxKefGcI/mYIeZ9WIg1bMLLekLsl+Bpih2i3wr7CTtwUWwTzWsZ\n3PE51FHsQIuH4XISuaGP4PorIq4mng3Da569PM1xyltiD28mhb5nm2hey1RLxg43ZYh/h+XeADuc\ne+OCKoIVu5JT6uqO0qt9pZS2WPjnuBJX9664kvhstCr2/MY98IHqfv1Xsdein9dj2QdjWAfiwmEc\nm9aOLbbcpJCe0dixM9F5bfoQYk9zZ1sbwupPkmeauSSwz0Lu3SjVLQuv6+bYledX3nXKx3Dpze22\nZ2eivD5rOjbkOpTXKq3vVPdMU+NHnDgQ4z4f6eru2FcY2ValQzVG61w85dmZKK9lldizh27Axlet\nR21P9Py+6gBYuJoflvbr7tgnVEMP1qFYJf+3Z3vReS17ErGjbtYbC4nO8nLvBe1X28Oy6/cEtths\nXYLYbm+iPpUaN89KmTd2Jsp7i01fOnx4e67YZUi1wFmm1HyJTnVf4sa7ka8xYuyZmqSs7lzVtfQ5\n1c2D/9Wdq5pXdWsaTroGLnW5Vs3fcL1zv7qYzzTxf9y2dOwos129U02iY7+VsEHwydJEZ/gMo2ue\nbaJ5LZN1k6AmMN8kz/Qp9cJ3GGNUCw/Eapr+hftUbu3p/I7szJNv6F9V1u+lw2OT+ldNnD7b96J5\nbVqNms4ztngvFMblCV/7OiaYbCQD5D1N/wX8wRK5TpgTOQAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$2.71828182845905^{- x}$$"
],
"text/plain": [
" -x\n",
"2.71828182845905 "
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAARBAMAAAARYzyGAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMlTvq5l2Zoki\nu0Rn3bgMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADUElEQVRIDcWVS0wTURSG/yktQ9spNLpA4qIl\nMTFoDES3Kk1ciHHRrsToAl9xYVDqE2iMVKKBBMNjJQQTutIgCxs1ASHGURNC3Ahx4UaSRpcueAho\nQBjPOXc6Q4C1nMXMud898//3zJ3eAtsQwd4Tt7bBFujRo0c3GWsdQzEHGu3TMjjWNDaWFar1jeSg\ncFl/p13oQsmI1gAV4ZoEdo3dBUTS0+N9aJcj5snkU7r7zsvAk9AOO7QUwTUevLAsKykVRcBtKJyG\nngBOhQEXSkYPHAEarBVoVzBuQiQLLSsL3zTFFFBiUo0dpW1zkj0FXuYZrgOXefAI0CEVZ4APCvsy\nMKJaUwsZu1AywH8V2N9sQj8HPQ2R9Aw/cWQRR8wd+JXxV6COtFRcBFpilGaB1yRHFW+BnRAcnIRW\nBQxQsQslAz4dAJL0XEkUwVmIpK4U+WrMvmLFfNjGv4GZRJ5VhpUxUNytjGemMArBxavhYEYZu1Ay\nWqhtHJmEfwEiuc5Ya/7ZH85bKFn6EP6QMZk4Ua1KQgR4abo1lOE5wi3Lo5Rwxy5UWcAk49TgAxRS\nx2tKUn/eG+MHN4fq2FgB4ll3tpjGHOwhFe/meSWMg1Y5ZWwMF0pWBjKewuccvWfPspIMhQ35Tql6\nQ9jG1LEYa2q6ICl3bZZuXOG7V32NUsa1P5ZyytiFKsuyMRCqQgXOLhp5yS8MndhxiGPPplcdiUpJ\nuyr084iNn8G/RD0SDiQxQL9A7tiFkvliyrjoLwJ9Jxec3esyldiGq+qYv4S6BE2pVv1pVVUwSXeu\nqKcXkgDjQhPeRWXsQsl2g4x13l56SJ9TkjdoQ3I03hy28TQdF9SEHadhmJxGknShCn7legaMI0RS\nYuzC9zJd29i41EzfVREtFJ4kRPIC0OXq0oQTtjH92lsd5k3Db/IonqaL3XHI9DIuNIFvTsc2pI5D\nxHGVT7OCaKAecZMPkFakgYM0sUWwLG1sQUK75Mx+HBu5yRB1aWJc8T2MQQguLkegWxm7UDIqXYVB\nJ5bpL9f20QJYcgK0iq2i6M58ijdWe/M45sxXWtYv2e3xBJ3IXOFtoj8JwTje0Qm0NezNrIMyTUeZ\nlcLE8H1gpD1HZwNLGh1NW79px+w/Jf8Ai2Q6WOYdK4QAAAAASUVORK5CYII=\n",
"text/latex": [
"$$- 1.0 \\cdot 2.71828182845905^{- x}$$"
],
"text/plain": [
" -x\n",
"-1.0⋅2.71828182845905 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sympy.diff(z)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
week13/.ipynb_checkpoints/Initial-Value-checkpoint.ipynb 0 → 100644
View file @ f7f9d76a
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Optimization\n",
"\n",
"> หาค่า $x$ ที่ทำให้ $F(x)$ มีค่าน้อยที่สุด และทำให้ $g(x) = 0, h(x) \\ge 0$\n",
"\n",
"$F(x)$ - objective function\n",
"\n",
"$x$ - design variable\n",
"\n",
"$g(x), h(x)$ - constraint functions\n",
"\n",
"minimize $F(x)$ = maximize $-F(x)$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> จงหาค่า $x$ ในช่วง $[a,b]$ ที่ทำ $f(x)$ มีค่าน้อยที่สุด\n",
"\n",
"## Bracketing\n",
"- $x = a$\n",
"- เลือกค่า $x$ ใหม่ ตามทิศทาง $dx = +\\Delta$ หรือ $dx = -\\Delta$ ที่ทำให้ค่า $f(x+dx) \\lt f(x)$\n",
"- กำหนดให้ $x = x + dx$ ไปเรื่อยๆ จนกว่า $f(x+dx) > f(x)$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Golden Section Search\n",
"\n",
"```python\n",
"def search(f,a,b,tol=1.0e-9):\n",
" nIter = int(math.ceil(-2.078087*math.log(tol/abs(b-a))))\n",
" R = 0.618033989\n",
" C = 1.0 - R\n",
" # First telescoping\n",
" x1 = R*a + C*b; x2 = C*a + R*b\n",
" f1 = f(x1); f2 = f(x2)\n",
" # Main loop\n",
" for i in range(nIter):\n",
" if f1 > f2:\n",
" a = x1\n",
" x1 = x2; f1 = f2\n",
" x2 = C*a + R*b; f2 = f(x2)\n",
" else:\n",
" b = x2\n",
" x2 = x1; f2 = f1\n",
" x1 = R*a + C*b; f1 = f(x1)\n",
"\n",
" if f1 < f2: return x1,f1\n",
" return x2,f2\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Powell's Method\n",
"\n",
"* Choose a point $x_0$ in the design space\n",
"* Choose the starting vectors $v_i$ , $i=1, 2, . . . , n$ (the usual choice is $v_i = e_i$ \n",
" * where $e_i$ is the unit vector in the $x_i$ -coordinate direction)\n",
"* Cycle:\n",
" * Loop over $i = 1, 2, . . . , n$:\n",
" * Minimize $F(x)$ along the line through $x_{i−1}$ in the direction of $v_i$\n",
" * Let the minimum point be $x_i$\n",
" * End loop.\n",
" * $v_{n+1} \\leftarrow x_0 − x_n$\n",
" * Minimize $F(x)$ along the line through $x_0$ in the direction of $v_{n+1}$\n",
" * Let the minimum point be $x_{n+1}$\n",
" * if $|x n+1 − x 0 | < ε$ exit loop\n",
" * Loop over $i = 1, 2, . . . , n$\n",
" * $v_i = v_{i+1}$ ($v_1$ is discarded; the other vectors are reused).\n",
" * End loop\n",
" * $x_0 = x_{n+1}$\n",
"* End cycle\n",
"\n",
"\n",
"```python\n",
"import numpy as np\n",
"import math\n",
"\n",
"def powell(F,x,h=0.1,tol=1.0e-6):\n",
"\n",
" def f(s): return F(x + s*v)\n",
"\n",
" n = len(x)\n",
" df = np.zeros(n)\n",
" u = np. identity(n)\n",
"\n",
" for j in range(30):\n",
" xOld = x.copy()\n",
" fOld = F(xOld)\n",
" \n",
" for i in range(n):\n",
" v = u[i]\n",
" a,b = bracket(f,0.0,h)\n",
" s,fMin = search(f,a,b)\n",
" df[i] = fOld - fMin\n",
" fOld = fMin\n",
" x = x + s*v\n",
"\n",
" v = x - xOld\n",
" a,b = bracket(f,0.0,h)\n",
" s,fLast = search(f,a,b)\n",
" x = x + s*v\n",
" if math.sqrt(np.dot(x-xOld,x-xOld)/n) < tol: return x,j+1\n",
" iMax = np.argmax(df)\n",
" for i in range(iMax,n-1):\n",
" u[i] = u[i+1]\n",
" u[n-1] = v\n",
"\n",
" print(\"Powell did not converge\")\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## example\n",
"1. จงเขียนโปรแกรมหา $x$ ที่ทำให้ $f(x) = 100(y-x^2)^2 + (1-x)^2$ ด้วยวิธี Powell โดยเริ่มจากจุด $(-1,1)$\n",
"\n",
"```python\n",
"from numpy import array\n",
"\n",
"def F(x): return 100.0*(x[1] - x[0]**2)**2 + (1 - x[0])**2\n",
"\n",
"xStart = array([-1.0, 1.0])\n",
"xMin,nIter = powell(F,xStart)\n",
"print(\"x =\",xMin)\n",
"print(\"F(x) =\",F(xMin))\n",
"print(\"Number of cycles =\",nIter)\n",
"input (\"Press return to exit\")\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Downhil Simplex Method (Nelder-Mead Method)\n",
"\n",
"* Choose a starting simplex\n",
"* Cycle until $d \\le \\epsilon$ \n",
" * Try reflection\n",
" * if new vertex $\\le$ old $Lo$: accept reflection\n",
" * Try expansion\n",
" * if new vertex $\\le$ old $Lo$: accept expansion\n",
" * else:\n",
" * if new vertex $>$ old $Hi$:\n",
" * Try contraction\n",
" * if new vertex $\\le$ old $Hi$: accept contraction\n",
" * else: use shrinkage\n",
" \n",
" \n",
"```python\n",
"import numpy as np\n",
"import math\n",
"\n",
"def downhill(F,xStart,side=0.1,tol=1.0e-6):\n",
" n = len(xStart)\n",
" x = np.zeros((n+1,n))\n",
" f = np.zeros(n+1)\n",
" x[0] = xStart\n",
" for i in range(1,n+1):\n",
" x[i] = xStart\n",
" x[i,i-1] = xStart[i-1] + side\n",
"\n",
" for i in range(n+1): f[i] = F(x[i])\n",
"\n",
" for k in range(500):\n",
" iLo = np.argmin(f)\n",
" iHi = np.argmax(f)\n",
" d = (-(n+1)*x[iHi] + np.sum(x,axis=0))/n\n",
" if math.sqrt(np.dot(d,d)/n) < tol: return x[iLo]\n",
"\n",
" # Try reflection\n",
" xNew = x[iHi] + 2.0*d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iLo]:\n",
" # Accept reflection\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
"\n",
" # Try expanding the reflection\n",
" xNew = x[iHi] + d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iLo]:\n",
" # Accept expansion\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" if fNew <= f[iHi]:\n",
" # Accept reflection\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" # Try contraction\n",
" xNew = x[iHi] + 0.5*d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iHi]: # Accept contraction\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" # Use shrinkage\n",
" for i in range(len(x)):\n",
" if i != iLo:395\n",
" x[i] = (x[i] - x[iLo])*0.5\n",
" f[i] = F(x[i])\n",
"\n",
" print(\"Too many iterations in downhill\")\n",
" return x[iLo]\n",
"```"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
week13/Initial-Value.ipynb 0 → 100644
View file @ f7f9d76a
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Optimization\n",
"\n",
"> หาค่า $x$ ที่ทำให้ $F(x)$ มีค่าน้อยที่สุด และทำให้ $g(x) = 0, h(x) \\ge 0$\n",
"\n",
"$F(x)$ - objective function\n",
"\n",
"$x$ - design variable\n",
"\n",
"$g(x), h(x)$ - constraint functions\n",
"\n",
"minimize $F(x)$ = maximize $-F(x)$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> จงหาค่า $x$ ในช่วง $[a,b]$ ที่ทำ $f(x)$ มีค่าน้อยที่สุด\n",
"\n",
"## Bracketing\n",
"- $x = a$\n",
"- เลือกค่า $x$ ใหม่ ตามทิศทาง $dx = +\\Delta$ หรือ $dx = -\\Delta$ ที่ทำให้ค่า $f(x+dx) \\lt f(x)$\n",
"- กำหนดให้ $x = x + dx$ ไปเรื่อยๆ จนกว่า $f(x+dx) > f(x)$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Golden Section Search\n",
"\n",
"```python\n",
"def search(f,a,b,tol=1.0e-9):\n",
" nIter = int(math.ceil(-2.078087*math.log(tol/abs(b-a))))\n",
" R = 0.618033989\n",
" C = 1.0 - R\n",
" # First telescoping\n",
" x1 = R*a + C*b; x2 = C*a + R*b\n",
" f1 = f(x1); f2 = f(x2)\n",
" # Main loop\n",
" for i in range(nIter):\n",
" if f1 > f2:\n",
" a = x1\n",
" x1 = x2; f1 = f2\n",
" x2 = C*a + R*b; f2 = f(x2)\n",
" else:\n",
" b = x2\n",
" x2 = x1; f2 = f1\n",
" x1 = R*a + C*b; f1 = f(x1)\n",
"\n",
" if f1 < f2: return x1,f1\n",
" return x2,f2\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Powell's Method\n",
"\n",
"* Choose a point $x_0$ in the design space\n",
"* Choose the starting vectors $v_i$ , $i=1, 2, . . . , n$ (the usual choice is $v_i = e_i$ \n",
" * where $e_i$ is the unit vector in the $x_i$ -coordinate direction)\n",
"* Cycle:\n",
" * Loop over $i = 1, 2, . . . , n$:\n",
" * Minimize $F(x)$ along the line through $x_{i−1}$ in the direction of $v_i$\n",
" * Let the minimum point be $x_i$\n",
" * End loop.\n",
" * $v_{n+1} \\leftarrow x_0 − x_n$\n",
" * Minimize $F(x)$ along the line through $x_0$ in the direction of $v_{n+1}$\n",
" * Let the minimum point be $x_{n+1}$\n",
" * if $|x n+1 − x 0 | < ε$ exit loop\n",
" * Loop over $i = 1, 2, . . . , n$\n",
" * $v_i = v_{i+1}$ ($v_1$ is discarded; the other vectors are reused).\n",
" * End loop\n",
" * $x_0 = x_{n+1}$\n",
"* End cycle\n",
"\n",
"\n",
"```python\n",
"import numpy as np\n",
"import math\n",
"\n",
"def powell(F,x,h=0.1,tol=1.0e-6):\n",
"\n",
" def f(s): return F(x + s*v)\n",
"\n",
" n = len(x)\n",
" df = np.zeros(n)\n",
" u = np. identity(n)\n",
"\n",
" for j in range(30):\n",
" xOld = x.copy()\n",
" fOld = F(xOld)\n",
" \n",
" for i in range(n):\n",
" v = u[i]\n",
" a,b = bracket(f,0.0,h)\n",
" s,fMin = search(f,a,b)\n",
" df[i] = fOld - fMin\n",
" fOld = fMin\n",
" x = x + s*v\n",
"\n",
" v = x - xOld\n",
" a,b = bracket(f,0.0,h)\n",
" s,fLast = search(f,a,b)\n",
" x = x + s*v\n",
" if math.sqrt(np.dot(x-xOld,x-xOld)/n) < tol: return x,j+1\n",
" iMax = np.argmax(df)\n",
" for i in range(iMax,n-1):\n",
" u[i] = u[i+1]\n",
" u[n-1] = v\n",
"\n",
" print(\"Powell did not converge\")\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## example\n",
"1. จงเขียนโปรแกรมหา $x$ ที่ทำให้ $f(x) = 100(y-x^2)^2 + (1-x)^2$ ด้วยวิธี Powell โดยเริ่มจากจุด $(-1,1)$\n",
"\n",
"```python\n",
"from numpy import array\n",
"\n",
"def F(x): return 100.0*(x[1] - x[0]**2)**2 + (1 - x[0])**2\n",
"\n",
"xStart = array([-1.0, 1.0])\n",
"xMin,nIter = powell(F,xStart)\n",
"print(\"x =\",xMin)\n",
"print(\"F(x) =\",F(xMin))\n",
"print(\"Number of cycles =\",nIter)\n",
"input (\"Press return to exit\")\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Downhil Simplex Method (Nelder-Mead Method)\n",
"\n",
"* Choose a starting simplex\n",
"* Cycle until $d \\le \\epsilon$ \n",
" * Try reflection\n",
" * if new vertex $\\le$ old $Lo$: accept reflection\n",
" * Try expansion\n",
" * if new vertex $\\le$ old $Lo$: accept expansion\n",
" * else:\n",
" * if new vertex $>$ old $Hi$:\n",
" * Try contraction\n",
" * if new vertex $\\le$ old $Hi$: accept contraction\n",
" * else: use shrinkage\n",
" \n",
" \n",
"```python\n",
"import numpy as np\n",
"import math\n",
"\n",
"def downhill(F,xStart,side=0.1,tol=1.0e-6):\n",
" n = len(xStart)\n",
" x = np.zeros((n+1,n))\n",
" f = np.zeros(n+1)\n",
" x[0] = xStart\n",
" for i in range(1,n+1):\n",
" x[i] = xStart\n",
" x[i,i-1] = xStart[i-1] + side\n",
"\n",
" for i in range(n+1): f[i] = F(x[i])\n",
"\n",
" for k in range(500):\n",
" iLo = np.argmin(f)\n",
" iHi = np.argmax(f)\n",
" d = (-(n+1)*x[iHi] + np.sum(x,axis=0))/n\n",
" if math.sqrt(np.dot(d,d)/n) < tol: return x[iLo]\n",
"\n",
" # Try reflection\n",
" xNew = x[iHi] + 2.0*d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iLo]:\n",
" # Accept reflection\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
"\n",
" # Try expanding the reflection\n",
" xNew = x[iHi] + d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iLo]:\n",
" # Accept expansion\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" if fNew <= f[iHi]:\n",
" # Accept reflection\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" # Try contraction\n",
" xNew = x[iHi] + 0.5*d\n",
" fNew = F(xNew)\n",
" if fNew <= f[iHi]: # Accept contraction\n",
" x[iHi] = xNew\n",
" f[iHi] = fNew\n",
" else:\n",
" # Use shrinkage\n",
" for i in range(len(x)):\n",
" if i != iLo:395\n",
" x[i] = (x[i] - x[iLo])*0.5\n",
" f[i] = F(x[i])\n",
"\n",
" print(\"Too many iterations in downhill\")\n",
" return x[iLo]\n",
"```"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment