Skip to content
Permalink
d7c884f285
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

ME3255F2017/06_roots and optimization/.ipynb_checkpoints/06_roots and optimization-checkpoint.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
913 lines (913 sloc) 43.4 KB
Raw Blame
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true,
"scrolled": true,
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"%plot --format svg"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Roots and Optimization\n",
"## Bracketing ch. 5"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Given a function, numerical or analytical, it's not always possible a given variable. \n",
"\n",
"### Freefall example:\n",
"If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n",
"\n",
"$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n",
"\n",
"Cannot solve for m \n",
"\n",
"Instead, solve the problem by creating a new function f(m) where\n",
"\n",
"$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n",
"\n",
"When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"scrolled": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"<svg height=\"420px\" viewBox=\"0 0 560 420\" width=\"560px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
"\n",
"<title>Gnuplot</title>\n",
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n",
"\n",
"<g id=\"gnuplot_canvas\">\n",
"\n",
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n",
"<defs>\n",
"\n",
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n",
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt6\" stroke=\"none\" xlink:href=\"#gpPt5\"/>\n",
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt8\" stroke=\"none\" xlink:href=\"#gpPt7\"/>\n",
"\t<use id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\" xlink:href=\"#gpPt7\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt10\" stroke=\"none\" xlink:href=\"#gpPt9\"/>\n",
"\t<use id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\" xlink:href=\"#gpPt3\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt12\" stroke=\"none\" xlink:href=\"#gpPt11\"/>\n",
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt14\" stroke=\"none\" xlink:href=\"#gpPt13\"/>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
"\t</filter>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
"\t</filter>\n",
"</defs>\n",
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n",
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"37.3,384.0 534.9,384.0 534.9,16.8 37.3,16.8 \"/>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,384.0 L49.8,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,390.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-5</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,322.8 L49.8,322.8 M535.0,322.8 L522.5,322.8 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,328.8)\">\n",
"\t\t<text><tspan font-family=\"{}\">-4</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,261.6 L49.8,261.6 M535.0,261.6 L522.5,261.6 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,267.6)\">\n",
"\t\t<text><tspan font-family=\"{}\">-3</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,200.3 L49.8,200.3 M535.0,200.3 L522.5,200.3 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,206.3)\">\n",
"\t\t<text><tspan font-family=\"{}\">-2</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,139.1 L49.8,139.1 M535.0,139.1 L522.5,139.1 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,145.1)\">\n",
"\t\t<text><tspan font-family=\"{}\">-1</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,77.9 L49.8,77.9 M535.0,77.9 L522.5,77.9 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,83.9)\">\n",
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,16.7 L49.8,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,22.7)\">\n",
"\t\t<text><tspan font-family=\"{}\">1</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M85.5,384.0 L85.5,371.5 M85.5,16.7 L85.5,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(85.5,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">60</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M149.7,384.0 L149.7,371.5 M149.7,16.7 L149.7,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(149.7,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">80</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M213.9,384.0 L213.9,371.5 M213.9,16.7 L213.9,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(213.9,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">100</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M278.1,384.0 L278.1,371.5 M278.1,16.7 L278.1,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(278.1,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">120</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M342.3,384.0 L342.3,371.5 M342.3,16.7 L342.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(342.3,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">140</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M406.6,384.0 L406.6,371.5 M406.6,16.7 L406.6,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(406.6,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">160</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M470.8,384.0 L470.8,371.5 M470.8,16.7 L470.8,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(470.8,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">180</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">200</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,16.7 L37.3,384.0 L535.0,384.0 L535.0,16.7 L37.3,16.7 Z \" stroke=\"black\"/></g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M53.4,358.3 L58.2,347.4 L63.1,337.1 L68.0,327.2 L72.8,317.8 L77.7,308.7 L82.5,300.0 L87.4,291.6 L92.3,283.6 L97.1,275.9 L102.0,268.4 L106.9,261.3 L111.7,254.4 L116.6,247.7 L121.5,241.3 L126.3,235.1 L131.2,229.1 L136.1,223.3 L140.9,217.7 L145.8,212.2 L150.7,207.0 L155.5,201.9 L160.4,197.0 L165.3,192.2 L170.1,187.5 L175.0,183.0 L179.8,178.6 L184.7,174.4 L189.6,170.3 L194.4,166.2 L199.3,162.3 L204.2,158.5 L209.0,154.8 L213.9,151.2 L218.8,147.7 L223.6,144.3 L228.5,140.9 L233.4,137.7 L238.2,134.5 L243.1,131.4 L248.0,128.4 L252.8,125.4 L257.7,122.6 L262.6,119.7 L267.4,117.0 L272.3,114.3 L277.1,111.7 L282.0,109.1 L286.9,106.6 L291.7,104.1 L296.6,101.7 L301.5,99.4 L306.3,97.1 L311.2,94.8 L316.1,92.6 L320.9,90.4 L325.8,88.3 L330.7,86.2 L335.5,84.2 L340.4,82.2 L345.3,80.2 L350.1,78.3 L355.0,76.4 L359.9,74.6 L364.7,72.8 L369.6,71.0 L374.5,69.2 L379.3,67.5 L384.2,65.8 L389.0,64.2 L393.9,62.5 L398.8,60.9 L403.6,59.4 L408.5,57.8 L413.4,56.3 L418.2,54.8 L423.1,53.3 L428.0,51.9 L432.8,50.5 L437.7,49.1 L442.6,47.7 L447.4,46.3 L452.3,45.0 L457.2,43.7 L462.0,42.4 L466.9,41.1 L471.8,39.9 L476.6,38.7 L481.5,37.5 L486.3,36.3 L491.2,35.1 L496.1,33.9 L500.9,32.8 L505.8,31.7 L510.7,30.6 L515.5,29.5 L520.4,28.4 L525.3,27.3 L530.1,26.3 L535.0,25.3 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n",
"\t</g>\n",
"\t<g id=\"gnuplot_plot_2a\"><title>gnuplot_plot_2a</title>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M53.4,77.9 L58.2,77.9 L63.1,77.9 L68.0,77.9 L72.8,77.9 L77.7,77.9 L82.5,77.9 L87.4,77.9 L92.3,77.9 L97.1,77.9 L102.0,77.9 L106.9,77.9 L111.7,77.9 L116.6,77.9 L121.5,77.9 L126.3,77.9 L131.2,77.9 L136.1,77.9 L140.9,77.9 L145.8,77.9 L150.7,77.9 L155.5,77.9 L160.4,77.9 L165.3,77.9 L170.1,77.9 L175.0,77.9 L179.8,77.9 L184.7,77.9 L189.6,77.9 L194.4,77.9 L199.3,77.9 L204.2,77.9 L209.0,77.9 L213.9,77.9 L218.8,77.9 L223.6,77.9 L228.5,77.9 L233.4,77.9 L238.2,77.9 L243.1,77.9 L248.0,77.9 L252.8,77.9 L257.7,77.9 L262.6,77.9 L267.4,77.9 L272.3,77.9 L277.1,77.9 L282.0,77.9 L286.9,77.9 L291.7,77.9 L296.6,77.9 L301.5,77.9 L306.3,77.9 L311.2,77.9 L316.1,77.9 L320.9,77.9 L325.8,77.9 L330.7,77.9 L335.5,77.9 L340.4,77.9 L345.3,77.9 L350.1,77.9 L355.0,77.9 L359.9,77.9 L364.7,77.9 L369.6,77.9 L374.5,77.9 L379.3,77.9 L384.2,77.9 L389.0,77.9 L393.9,77.9 L398.8,77.9 L403.6,77.9 L408.5,77.9 L413.4,77.9 L418.2,77.9 L423.1,77.9 L428.0,77.9 L432.8,77.9 L437.7,77.9 L442.6,77.9 L447.4,77.9 L452.3,77.9 L457.2,77.9 L462.0,77.9 L466.9,77.9 L471.8,77.9 L476.6,77.9 L481.5,77.9 L486.3,77.9 L491.2,77.9 L496.1,77.9 L500.9,77.9 L505.8,77.9 L510.7,77.9 L515.5,77.9 L520.4,77.9 L525.3,77.9 L530.1,77.9 L535.0,77.9 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n",
"\t</g>\n",
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 128, 0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"</g>\n",
"</svg>"
],
"text/plain": [
"<IPython.core.display.SVG object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"setdefaults\n",
"g=9.81; % acceleration due to gravity\n",
"m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n",
"c_d=0.25; % drag coefficient\n",
"t=4; % at time = 4 seconds\n",
"v=36; % speed must be 36 m/s\n",
"f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n",
"\n",
"plot(m,f_m(m),m,zeros(length(m),1))\n",
"axis([45 200 -5 1])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"scrolled": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans = -0.056985\r\n"
]
}
],
"source": [
"f_m(140)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n",
"\n",
"Better methods are the \n",
"1. Bracketing methods\n",
"2. Open methods\n",
"\n",
"Both need an initial guess. \n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Incremental method (Brute force)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"^C\r\n"
]
}
],
"source": [
"function xb = incsearch(func,xmin,xmax,ns)\n",
"% incsearch: incremental search root locator\n",
"% xb = incsearch(func,xmin,xmax,ns):\n",
"% finds brackets of x that contain sign changes\n",
"% of a function on an interval\n",
"% input:\n",
"% func = name of function\n",
"% xmin, xmax = endpoints of interval\n",
"% ns = number of subintervals (default = 50)\n",
"% output:\n",
"% xb(k,1) is the lower bound of the kth sign change\n",
"% xb(k,2) is the upper bound of the kth sign change\n",
"% If no brackets found, xb = [].\n",
"if nargin < 3, error('at least 3 arguments required'), end\n",
"if nargin < 4, ns = 50; end %if ns blank set to 50\n",
"% Incremental search\n",
"x = linspace(xmin,xmax,ns);\n",
"f = func(x);\n",
"nb = 0; xb = []; %xb is null unless sign change detected\n",
"for k = 1:length(x)-1\n",
" if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n",
" nb = nb + 1;\n",
" xb(nb,1) = x(k);\n",
" xb(nb,2) = x(k+1);\n",
" end\n",
"end\n",
"if isempty(xb) %display that no brackets were found\n",
" fprintf('no brackets found\\n')\n",
" fprintf('check interval or increase ns\\n')\n",
"else\n",
" fprintf('number of brackets: %i\\n',nb) %display number of brackets\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"'incsearch' is a function from the file /home/ryan/Documents/UConn/ME3255/me3255_S2017/course_git/lecture_06/incsearch.m\n",
"\n",
" incsearch: incremental search root locator\n",
" xb = incsearch(func,xmin,xmax,ns):\n",
" finds brackets of x that contain sign changes\n",
" of a function on an interval\n",
" input:\n",
" func = name of function\n",
" xmin, xmax = endpoints of interval\n",
" ns = number of subintervals (default = 50)\n",
" output:\n",
" xb(k,1) is the lower bound of the kth sign change\n",
" xb(k,2) is the upper bound of the kth sign change\n",
" If no brackets found, xb = [].\n",
"\n",
"\n",
"Additional help for built-in functions and operators is\n",
"available in the online version of the manual. Use the command\n",
"'doc <topic>' to search the manual index.\n",
"\n",
"Help and information about Octave is also available on the WWW\n",
"at http://www.octave.org and via the help@octave.org\n",
"mailing list.\n"
]
}
],
"source": [
"help incsearch"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false,
"scrolled": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of brackets: 1\n",
"ans =\n",
"\n",
" 142.42 143.94\n",
"\n"
]
}
],
"source": [
"incsearch(f_m,50, 200,100)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Bisection method"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Divide interval in half until error is reduced to some level\n",
"\n",
"in previous example of freefall, choose x_l=50, x_u=200\n",
"\n",
"x_r = (50+200)/2 = 125\n",
"\n",
"f_m(125) = -0.408\n",
"\n",
"x_r= (125+200)/2 = 162.5\n",
"\n",
"f_m(162.5) = 0.3594\n",
"\n",
"x_r = (125+162.5)/2=143.75\n",
"\n",
"f_m(143.75)= 0.0206"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"scrolled": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans = 0.020577\r\n"
]
}
],
"source": [
"f_m(143.75)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Bisect"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Much better root locator, with 4 iterations, our function is already close to zero\n",
"\n",
"Automate this with a function:\n",
"`bisect.m`"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"^C\r\n"
]
}
],
"source": [
"function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n",
"% bisect: root location zeroes\n",
"% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
"% uses bisection method to find the root of func\n",
"% input:\n",
"% func = name of function\n",
"% xl, xu = lower and upper guesses\n",
"% es = desired relative error (default = 0.0001%)\n",
"% maxit = maximum allowable iterations (default = 50)\n",
"% p1,p2,... = additional parameters used by func\n",
"% output:\n",
"% root = real root\n",
"% fx = function value at root\n",
"% ea = approximate relative error (%)\n",
"% iter = number of iterations\n",
"if nargin<3,error('at least 3 input arguments required'),end\n",
"test = func(xl,varargin{:})*func(xu,varargin{:});\n",
"if test>0,error('no sign change'),end\n",
"if nargin<4|isempty(es), es=0.0001;end\n",
"if nargin<5|isempty(maxit), maxit=50;end\n",
"iter = 0; xr = xl; ea = 100;\n",
"while (1)\n",
" xrold = xr;\n",
" xr = (xl + xu)/2;\n",
" iter = iter + 1;\n",
" if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n",
" test = func(xl,varargin{:})*func(xr,varargin{:});\n",
" if test < 0\n",
" xu = xr;\n",
" elseif test > 0\n",
" xl = xr;\n",
" else\n",
" ea = 0;\n",
" end\n",
" if ea <= es | iter >= maxit,break,end\n",
"end\n",
"root = xr; fx = func(xr, varargin{:});"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"'bisect' is a function from the file /home/ryan/Documents/UConn/ME3255/course_git_F2017/06_roots and optimization/bisect.m\n",
"\n",
" bisect: root location zeroes\n",
" [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
" uses bisection method to find the root of func\n",
" input:\n",
" func = name of function\n",
" xl, xu = lower and upper guesses\n",
" es = desired relative error (default = 0.0001%)\n",
" maxit = maximum allowable iterations (default = 50)\n",
" p1,p2,... = additional parameters used by func\n",
" output:\n",
" root = real root\n",
" fx = function value at root\n",
" ea = approximate relative error (%)\n",
" iter = number of iterations\n",
"\n",
"\n",
"Additional help for built-in functions and operators is\n",
"available in the online version of the manual. Use the command\n",
"'doc <topic>' to search the manual index.\n",
"\n",
"Help and information about Octave is also available on the WWW\n",
"at http://www.octave.org and via the help@octave.org\n",
"mailing list.\n"
]
}
],
"source": [
"help bisect"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## False position (linear interpolation)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n",
"\n",
"$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false,
"scrolled": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"<svg height=\"420px\" viewBox=\"0 0 560 420\" width=\"560px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
"\n",
"<title>Gnuplot</title>\n",
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n",
"\n",
"<g id=\"gnuplot_canvas\">\n",
"\n",
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n",
"<defs>\n",
"\n",
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n",
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt6\" stroke=\"none\" xlink:href=\"#gpPt5\"/>\n",
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt8\" stroke=\"none\" xlink:href=\"#gpPt7\"/>\n",
"\t<use id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\" xlink:href=\"#gpPt7\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt10\" stroke=\"none\" xlink:href=\"#gpPt9\"/>\n",
"\t<use id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\" xlink:href=\"#gpPt3\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt12\" stroke=\"none\" xlink:href=\"#gpPt11\"/>\n",
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt14\" stroke=\"none\" xlink:href=\"#gpPt13\"/>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
"\t</filter>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
"\t</filter>\n",
"</defs>\n",
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n",
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"37.3,384.0 534.9,384.0 534.9,16.8 37.3,16.8 \"/>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,384.0 L49.8,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,390.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-5</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,322.8 L49.8,322.8 M535.0,322.8 L522.5,322.8 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,328.8)\">\n",
"\t\t<text><tspan font-family=\"{}\">-4</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,261.6 L49.8,261.6 M535.0,261.6 L522.5,261.6 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,267.6)\">\n",
"\t\t<text><tspan font-family=\"{}\">-3</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,200.3 L49.8,200.3 M535.0,200.3 L522.5,200.3 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,206.3)\">\n",
"\t\t<text><tspan font-family=\"{}\">-2</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,139.1 L49.8,139.1 M535.0,139.1 L522.5,139.1 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,145.1)\">\n",
"\t\t<text><tspan font-family=\"{}\">-1</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,77.9 L49.8,77.9 M535.0,77.9 L522.5,77.9 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,83.9)\">\n",
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,16.7 L49.8,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(29.0,22.7)\">\n",
"\t\t<text><tspan font-family=\"{}\">1</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,384.0 L37.3,371.5 M37.3,16.7 L37.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(37.3,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M161.7,384.0 L161.7,371.5 M161.7,16.7 L161.7,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(161.7,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">50</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M286.2,384.0 L286.2,371.5 M286.2,16.7 L286.2,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(286.2,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">100</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M410.6,384.0 L410.6,371.5 M410.6,16.7 L410.6,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(410.6,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">150</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">200</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M37.3,16.7 L37.3,384.0 L535.0,384.0 L535.0,16.7 L37.3,16.7 Z \" stroke=\"black\"/></g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M161.7,358.3 L165.5,347.4 L169.3,337.1 L173.0,327.2 L176.8,317.8 L180.6,308.7 L184.3,300.0 L188.1,291.6 L191.9,283.6 L195.7,275.9 L199.4,268.4 L203.2,261.3 L207.0,254.4 L210.7,247.7 L214.5,241.3 L218.3,235.1 L222.1,229.1 L225.8,223.3 L229.6,217.7 L233.4,212.2 L237.1,207.0 L240.9,201.9 L244.7,197.0 L248.4,192.2 L252.2,187.5 L256.0,183.0 L259.8,178.6 L263.5,174.4 L267.3,170.3 L271.1,166.2 L274.8,162.3 L278.6,158.5 L282.4,154.8 L286.2,151.2 L289.9,147.7 L293.7,144.3 L297.5,140.9 L301.2,137.7 L305.0,134.5 L308.8,131.4 L312.5,128.4 L316.3,125.4 L320.1,122.6 L323.9,119.7 L327.6,117.0 L331.4,114.3 L335.2,111.7 L338.9,109.1 L342.7,106.6 L346.5,104.1 L350.2,101.7 L354.0,99.4 L357.8,97.1 L361.6,94.8 L365.3,92.6 L369.1,90.4 L372.9,88.3 L376.6,86.2 L380.4,84.2 L384.2,82.2 L388.0,80.2 L391.7,78.3 L395.5,76.4 L399.3,74.6 L403.0,72.8 L406.8,71.0 L410.6,69.2 L414.3,67.5 L418.1,65.8 L421.9,64.2 L425.7,62.5 L429.4,60.9 L433.2,59.4 L437.0,57.8 L440.7,56.3 L444.5,54.8 L448.3,53.3 L452.1,51.9 L455.8,50.5 L459.6,49.1 L463.4,47.7 L467.1,46.3 L470.9,45.0 L474.7,43.7 L478.4,42.4 L482.2,41.1 L486.0,39.9 L489.8,38.7 L493.5,37.5 L497.3,36.3 L501.1,35.1 L504.8,33.9 L508.6,32.8 L512.4,31.7 L516.1,30.6 L519.9,29.5 L523.7,28.4 L527.5,27.3 L531.2,26.3 L535.0,25.3 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n",
"\t</g>\n",
"\t<g id=\"gnuplot_plot_2a\"><title>gnuplot_plot_2a</title>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<g onmousemove=\"gnuplot_svg.showHypertext(evt,'')\" onmouseout=\"gnuplot_svg.hideHypertext()\"><title> </title>\n",
"\t<use color=\"rgb( 0, 128, 0)\" transform=\"translate(161.7,358.3) scale(12.00)\" xlink:href=\"#gpPt3\"/></g>\n",
"</g>\n",
"\t</g>\n",
"\t<g id=\"gnuplot_plot_3a\"><title>gnuplot_plot_3a</title>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<g onmousemove=\"gnuplot_svg.showHypertext(evt,'')\" onmouseout=\"gnuplot_svg.hideHypertext()\"><title> </title>\n",
"\t<use color=\"rgb(255, 0, 0)\" transform=\"translate(402.3,73.1) scale(12.00)\" xlink:href=\"#gpPt3\"/></g>\n",
"</g>\n",
"\t</g>\n",
"\t<g id=\"gnuplot_plot_4a\"><title>gnuplot_plot_4a</title>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M398.3,77.9 \" stroke=\"rgb( 0, 191, 191)\"/>\t<g onmousemove=\"gnuplot_svg.showHypertext(evt,'')\" onmouseout=\"gnuplot_svg.hideHypertext()\"><title> </title>\n",
"\t<use color=\"rgb( 0, 191, 191)\" transform=\"translate(398.3,77.9) scale(4.00)\" xlink:href=\"#gpPt6\"/></g>\n",
"</g>\n",
"\t</g>\n",
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 191, 191)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"</g>\n",
"</svg>"
],
"text/plain": [
"<IPython.core.display.SVG object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%xl=50; xu=200; \n",
"xu=xr;\n",
"xl=xl;\n",
"xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n",
"plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## False Position"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Much better root locator, with 4 iterations, our function is already close to zero\n",
"\n",
"Automate this with a function:\n",
"`falsepos.m`"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"^C\r\n"
]
}
],
"source": [
"function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n",
"% falsepos: root location zeroes\n",
"% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
"% uses false position method to find the root of func\n",
"% input:\n",
"% func = name of function\n",
"% xl, xu = lower and upper guesses\n",
"% es = desired relative error (default = 0.0001%)\n",
"% maxit = maximum allowable iterations (default = 50)\n",
"% p1,p2,... = additional parameters used by func\n",
"% output:\n",
"% root = real root\n",
"% fx = function value at root\n",
"% ea = approximate relative error (%)\n",
"% iter = number of iterations\n",
"if nargin<3,error('at least 3 input arguments required'),end\n",
"test = func(xl,varargin{:})*func(xu,varargin{:});\n",
"if test>0,error('no sign change'),end\n",
"if nargin<4|isempty(es), es=0.0001;end\n",
"if nargin<5|isempty(maxit), maxit=50;end\n",
"iter = 0; xr = xl; ea = 100;\n",
"while (1)\n",
" xrold = xr;\n",
" % xr = (xl + xu)/2; % bisect method\n",
" xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n",
" iter = iter + 1;\n",
" if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n",
" test = func(xl,varargin{:})*func(xr,varargin{:});\n",
" if test < 0\n",
" xu = xr;\n",
" elseif test > 0\n",
" xl = xr;\n",
" else\n",
" ea = 0;\n",
" end\n",
" if ea <= es | iter >= maxit,break,end\n",
"end\n",
"root = xr; fx = func(xr, varargin{:});"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "0.19.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}