{"id":188,"date":"2009-06-16T18:02:41","date_gmt":"2009-06-16T18:02:41","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/16\/rooting-around-in-matlab-part-2\/"},"modified":"2009-06-15T15:08:34","modified_gmt":"2009-06-15T15:08:34","slug":"rooting-around-in-matlab-part-2","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/16\/rooting-around-in-matlab-part-2\/","title":{"rendered":"Rooting Around in MATLAB &#8211; Part 2"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>In a <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=187\">recent post<\/a>, I started a series of posts about finding roots of equations, and how you might use MATLAB to teach some of this topic.\r\n          In particular, I will introduce the idea of <a href=\"http:\/\/en.wikipedia.org\/wiki\/Fixed_point_iteration\">fixed point iteration<\/a>.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Set Up<\/a><\/li>\r\n         <li><a href=\"#2\">Example<\/a><\/li>\r\n         <li><a href=\"#3\">Fixed Point Algorithm<\/a><\/li>\r\n         <li><a href=\"#8\">Let's Try to \"Zero\" in on the Solution<\/a><\/li>\r\n         <li><a href=\"#10\">First Iteration<\/a><\/li>\r\n         <li><a href=\"#11\">Second Iteration<\/a><\/li>\r\n         <li><a href=\"#12\">Third Iteration<\/a><\/li>\r\n         <li><a href=\"#13\">The Conundrum<\/a><\/li>\r\n         <li><a href=\"#14\">Confession<\/a><\/li>\r\n         <li><a href=\"#15\">The Series of Posts<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Set Up<a name=\"1\"><\/a><\/h3><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">clear <span style=\"color: #A020F0\">all<\/span>\r\nclose <span style=\"color: #A020F0\">all<\/span><\/pre><h3>Example<a name=\"2\"><\/a><\/h3>\r\n   <p>Let me restate the example function. We'll start with a simple cubic polynomial <tt>f<\/tt><\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_eq61733.png\"> <\/p>\r\n   <p>which I define as an anonymous function.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">f = @(x) x.^3 + x - 1;\r\nfplot(f,[-2 2])\r\ntitle <span style=\"color: #A020F0\">f<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_01.png\"> <h3>Fixed Point Algorithm<a name=\"3\"><\/a><\/h3>\r\n   <p>One way to find a solution to <tt>f(x)<\/tt> = 0 is to define another function, <tt>g<\/tt> such that <tt>g(x) = x<\/tt> where the function <tt>g<\/tt> is related to <tt>f<\/tt>.  If we can find a value <tt>x<\/tt> that solves the equation with <tt>g<\/tt>, this <b>fixed point<\/b> is also the zero of the function <tt>f<\/tt>.  How do we define such a suitable <tt>g<\/tt>? Looking at the equation for <tt>f<\/tt> and setting its value to zero, we can derive an equation for <tt>x<\/tt>.\r\n   <\/p>\r\n   <p>Here's one way we can rewrite <tt>f(x) = 0<\/tt> so that <tt>x<\/tt> is alone on one side of the equation.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_eq35464.png\"> <\/p>\r\n   <p>Let's call this function <tt>g1<\/tt>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">g1 = @(x) 1 - x.^3;\r\nfplot(g1,[-2 2]);\r\ng1str = <span style=\"color: #A020F0\">'1-x^3'<\/span>;\r\ntitle <span style=\"color: #A020F0\">g1<\/span>\r\ngrid <span style=\"color: #A020F0\">on<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_02.png\"> <p>Now let's add a straight line to represent<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_eq90685.png\"> <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">straightLine = @(x) x;\r\ntitle <span style=\"color: #A020F0\">''<\/span>\r\nhold <span style=\"color: #A020F0\">on<\/span>\r\nfplot(straightLine, [-2 2], <span style=\"color: #A020F0\">'g'<\/span>)\r\nlegend(g1str,<span style=\"color: #A020F0\">'x'<\/span>,<span style=\"color: #A020F0\">'Location'<\/span>,<span style=\"color: #A020F0\">'SouthEast'<\/span>)\r\n\r\naxis([-2 2 -2 2]), axis <span style=\"color: #A020F0\">equal<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_03.png\"> <p>The intersection of these 2 curves is a fixed point of <tt>g1(x)<\/tt>, and therefore a zero of <tt>f(x)<\/tt>.  Let's zoom in a bit.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">axis([0 1 0 1])<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_04.png\"> <p>That intersection looks like it happens around the same location as the answers we got from using <tt>roots<\/tt> and <tt>fzero<\/tt> in the previous post.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">fzsolution = fzero(f,0.5)<\/pre><pre style=\"font-style:oblique\">fzsolution =\r\n      0.68233\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">plot(fzsolution, g1(fzsolution),<span style=\"color: #A020F0\">'b*'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_05.png\"> <h3>Let's Try to \"Zero\" in on the Solution<a name=\"8\"><\/a><\/h3>\r\n   <p>First, we'll select the same starting value that we used when we called <tt>fzero<\/tt>.  The starting point is then (0.5,0.5).  And then we can easily calculate the next value of <tt>g1<\/tt>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">x(1) = 0.5;\r\ny(1) = x(1);\r\nx(2) = x(1);\r\ny(2) = g1(x(2));<\/pre><p>Let's plot the trajectory of the solutions, noting the first point with a red star.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">plot(x,y,<span style=\"color: #A020F0\">'r'<\/span>,x(1),y(1),<span style=\"color: #A020F0\">'r*'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_06.png\"> <h3>First Iteration<a name=\"10\"><\/a><\/h3>\r\n   <p>Now let's get the next estimate.  We'll set the current value of <tt>g1<\/tt> to the new <tt>x<\/tt> value, sliding onto the straight line, and then calculate the next value of <tt>g1<\/tt> from this new value for <tt>x<\/tt>.  And plot it.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">x(3) = y(2);\r\ny(3) = y(2);\r\nx(4) = x(3);\r\ny(4) = g1(x(4));\r\nplot(x,y,<span style=\"color: #A020F0\">'r'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_07.png\"> <h3>Second Iteration<a name=\"11\"><\/a><\/h3>\r\n   <p>Let's do the same set of steps again.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">n = 5;\r\nx(n) = y(n-1);\r\ny(n) = y(n-1);\r\nx(n+1) = x(n);\r\ny(n+1) = g1(x(n+1));\r\nplot(x,y,<span style=\"color: #A020F0\">'r'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_08.png\"> <h3>Third Iteration<a name=\"12\"><\/a><\/h3>\r\n   <p>And again.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">n = 7;\r\nx(n) = y(n-1);\r\ny(n) = y(n-1);\r\nx(n+1) = x(n);\r\ny(n+1) = g1(x(n+1));\r\nplot(x,y,<span style=\"color: #A020F0\">'r'<\/span>)\r\nhold <span style=\"color: #A020F0\">off<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/188\/findFP2_09.png\"> <h3>The Conundrum<a name=\"13\"><\/a><\/h3>\r\n   <p>By now, you have surely noticed that instead of spiraling into the solution we found with <tt>fzero<\/tt>, our estimates are spiraling out!  In the next post, we'll explore another definition for <tt>g(x) = x<\/tt> and see what happens from there.\r\n   <\/p>\r\n   <h3>Confession<a name=\"14\"><\/a><\/h3>\r\n   <p>Yes, I know in this post I kept overplotting lines successively.  But I didn't want us to get distracted by managing the plots.<\/p>\r\n   <h3>The Series of Posts<a name=\"15\"><\/a><\/h3>\r\n   <p>In addition to this post and the previous one, there will be one more. Until it is published, the 3rd link will not be useful.<\/p>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=187\">Rooting Around in MATLAB - Part 1<\/a><\/li>\r\n         <li><a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=188\">Rooting Around in MATLAB - Part 2 (this post)<\/a><\/li>\r\n         <li><a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=189\">Rooting Around in MATLAB - Part 3<\/a><\/li>\r\n      <\/ul>\r\n   <\/div><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_0b426223f51148dd8b3d38102cc95242() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='0b426223f51148dd8b3d38102cc95242 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 0b426223f51148dd8b3d38102cc95242';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2009 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_0b426223f51148dd8b3d38102cc95242()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.8<br><\/p>\r\n<\/div>\r\n<!--\r\n0b426223f51148dd8b3d38102cc95242 ##### SOURCE BEGIN #####\r\n%% Rooting Around in MATLAB - Part 2\r\n% In a <https:\/\/blogs.mathworks.com\/loren\/?p=187 recent post>, I started a\r\n% series of posts about finding roots of equations, and how you might use\r\n% MATLAB to teach some of this topic.  In particular, I will introduce the\r\n% idea of \r\n% <http:\/\/en.wikipedia.org\/wiki\/Fixed_point_iteration fixed point iteration>.\r\n%% Set Up\r\nclear all\r\nclose all\r\n%% Example\r\n% Let me restate the example function. We'll start with a simple cubic\r\n% polynomial |f|\r\n%\r\n% $f(x) = x^3 + x - 1$ \r\n%\r\n% which I define as an anonymous function.\r\nf = @(x) x.^3 + x - 1;\r\nfplot(f,[-2 2])\r\ntitle f\r\ngrid on\r\n%% Fixed Point Algorithm\r\n% One way to find a solution to |f(x)| = 0 is to define another function,\r\n% |g| such that |g(x) = x| where the function |g| is related to |f|.  If we\r\n% can find a value |x| that solves the equation with |g|, this *fixed\r\n% point* is also the zero of the function |f|.  How do we define such a\r\n% suitable |g|? Looking at the equation for |f| and setting its value to\r\n% zero, we can derive an equation for |x|.\r\n%\r\n% Here's one way we can rewrite |f(x) = 0| so that |x| is alone on one side\r\n% of the equation.\r\n%\r\n% $x = 1 - x^3$\r\n%\r\n% Let's call this function |g1|. \r\ng1 = @(x) 1 - x.^3;\r\nfplot(g1,[-2 2]);\r\ng1str = '1-x^3';\r\ntitle g1\r\ngrid on\r\n%%\r\n% Now let's add a straight line to represent \r\n%\r\n% $y = x$\r\n%\r\nstraightLine = @(x) x;\r\ntitle ''\r\nhold on\r\nfplot(straightLine, [-2 2], 'g')\r\nlegend(g1str,'x','Location','SouthEast')\r\n\r\naxis([-2 2 -2 2]), axis equal\r\n%%\r\n% The intersection of these 2 curves is a fixed point of |g1(x)|, and\r\n% therefore a zero of |f(x)|.  Let's zoom in a bit.\r\naxis([0 1 0 1])\r\n%%\r\n% That intersection looks like it happens around the same location as the\r\n% answers we got from using |roots| and |fzero| in the previous post.\r\nfzsolution = fzero(f,0.5)\r\n%%\r\nplot(fzsolution, g1(fzsolution),'b*')\r\n%% Let's Try to \"Zero\" in on the Solution\r\n% First, we'll select the same starting value that we used when we called\r\n% |fzero|.  The starting point is then (0.5,0.5).  And then we can easily\r\n% calculate the next value of |g1|.\r\nx(1) = 0.5;\r\ny(1) = x(1);\r\nx(2) = x(1);\r\ny(2) = g1(x(2));\r\n%%\r\n% Let's plot the trajectory of the solutions, noting the first\r\n% point with a red star.\r\nplot(x,y,'r',x(1),y(1),'r*')\r\n%% First Iteration\r\n% Now let's get the next estimate.  We'll set the current value of |g1| to\r\n% the new |x| value, sliding onto the straight line, and then calculate the\r\n% next value of |g1| from this new value for |x|.  And plot it.\r\nx(3) = y(2);\r\ny(3) = y(2);\r\nx(4) = x(3);\r\ny(4) = g1(x(4));\r\nplot(x,y,'r')\r\n%% Second Iteration\r\n% Let's do the same set of steps again.\r\nn = 5;\r\nx(n) = y(n-1);\r\ny(n) = y(n-1);\r\nx(n+1) = x(n);\r\ny(n+1) = g1(x(n+1));\r\nplot(x,y,'r')\r\n%% Third Iteration\r\n% And again.\r\nn = 7;\r\nx(n) = y(n-1);\r\ny(n) = y(n-1);\r\nx(n+1) = x(n);\r\ny(n+1) = g1(x(n+1));\r\nplot(x,y,'r')\r\nhold off\r\n%% The Conundrum\r\n% By now, you have surely noticed that instead of spiraling into the\r\n% solution we found with |fzero|, our estimates are spiraling out!  In the\r\n% next post, we'll explore another definition for |g(x) = x| and see what\r\n% happens from there.\r\n%% Confession\r\n% Yes, I know in this post I kept overplotting lines successively.  But I\r\n% didn't want us to get distracted by managing the plots.\r\n%% The Series of Posts\r\n% In addition to this post and the previous one, there will be one more.\r\n% Until it is published, the 3rd link will not be useful.\r\n%\r\n% * <https:\/\/blogs.mathworks.com\/loren\/?p=187 Rooting Around in MATLAB - Part 1>\r\n% * <https:\/\/blogs.mathworks.com\/loren\/?p=188 Rooting Around in MATLAB - Part 2 (this post)>\r\n% * <https:\/\/blogs.mathworks.com\/loren\/?p=189 Rooting Around in MATLAB - Part 3>\r\n%\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n##### SOURCE END ##### 0b426223f51148dd8b3d38102cc95242\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      In a recent post, I started a series of posts about finding roots of equations, and how you might use MATLAB to teach some of this topic.\r\n          In particular, I will... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/16\/rooting-around-in-matlab-part-2\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[25,31],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/188"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=188"}],"version-history":[{"count":0,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/188\/revisions"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=188"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=188"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=188"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}