{"id":189,"date":"2009-06-23T18:08:34","date_gmt":"2009-06-23T18:08:34","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/23\/rooting-around-in-matlab-part-3\/"},"modified":"2009-06-15T15:08:49","modified_gmt":"2009-06-15T15:08:49","slug":"rooting-around-in-matlab-part-3","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/23\/rooting-around-in-matlab-part-3\/","title":{"rendered":"Rooting Around in MATLAB &#8211; Part 3"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>I recently wrote a pair of posts (<a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=187\">1<\/a> and <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=188\">2<\/a>) about finding roots of equations, and how you might use MATLAB to explore this topic with emphasis on the method 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 Reminder<\/a><\/li>\r\n         <li><a href=\"#3\">First Fixed Point Iteration Attempt<\/a><\/li>\r\n         <li><a href=\"#5\">Second Way to Rewrite Equation<\/a><\/li>\r\n         <li><a href=\"#6\">Let's Iterate<\/a><\/li>\r\n         <li><a href=\"#7\">Let's Try to \"Zero\" in on the Solution<\/a><\/li>\r\n         <li><a href=\"#9\">First Iteration<\/a><\/li>\r\n         <li><a href=\"#10\">Second Iteration<\/a><\/li>\r\n         <li><a href=\"#11\">Third Iteration<\/a><\/li>\r\n         <li><a href=\"#12\">A Few More Rounds<\/a><\/li>\r\n         <li><a href=\"#13\">Fixed Point Found!<\/a><\/li>\r\n         <li><a href=\"#14\">Wrapping Up 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 Reminder<a name=\"2\"><\/a><\/h3>\r\n   <p>Let me restate the example function. Let's start with a simple cubic polynomial <tt>f(x)<\/tt><\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/189\/findFP3_eq61733.png\"> <\/p>\r\n   <p>which we define using 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\/189\/findFP3_01.png\"> <h3>First Fixed Point Iteration Attempt<a name=\"3\"><\/a><\/h3>\r\n   <p>Previously we rewrote <tt>f(x)=0<\/tt> so that <tt>x<\/tt> was 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\/189\/findFP3_eq35464.png\"> <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">g1 = @(x) 1 - x.^3;<\/pre><p>This function, <tt>g1(x)<\/tt> did not help find the root between 0 and 1 - every step took us further away from the solutions we found with <tt>roots<\/tt> and <tt>fzero<\/tt>.\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><h3>Second Way to Rewrite Equation<a name=\"5\"><\/a><\/h3>\r\n   <p>There's another way we can write the equation for the fixed point. Remember, we want to rearrange <tt>f(x)=0<\/tt> into something like <tt>g2(x)=x<\/tt>. We have already tried <tt>g1 = 1 - x^3<\/tt>.  We can also rewrite the equation as\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/189\/findFP3_eq32690.png\"> <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">g2 = @(x) 1.\/(x.^2+1);\r\nfplot(g2,[0 1]);\r\nhold <span style=\"color: #A020F0\">on<\/span>\r\nstraightLine = @(x) x;\r\nfplot(straightLine, [0 1], <span style=\"color: #A020F0\">'g'<\/span>)\r\nlegend(<span style=\"color: #A020F0\">'g2'<\/span>,<span style=\"color: #A020F0\">'x'<\/span>,<span style=\"color: #A020F0\">'Location'<\/span>,<span style=\"color: #A020F0\">'SouthEast'<\/span>)\r\ngrid <span style=\"color: #A020F0\">on<\/span>\r\naxis <span style=\"color: #A020F0\">equal<\/span>, axis([0 1 0 1])<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/189\/findFP3_02.png\"> <h3>Let's Iterate<a name=\"6\"><\/a><\/h3>\r\n   <p>Following the same idea from the last post, we now iterate, finding successive guesses of <tt>x<\/tt>, computing the value of <tt>g2(x)<\/tt>, setting this value to be the next estimate of <tt>x<\/tt>, etc.\r\n   <\/p>\r\n   <h3>Let's Try to \"Zero\" in on the Solution<a name=\"7\"><\/a><\/h3>\r\n   <p>First, we'll select a starting value <tt>x = 0.5<\/tt>, and <tt>y<\/tt>, from the straight line also at 0.5.\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) = g2(x(2));<\/pre><p>And let's plot the trajectory of the solutions, starting with the first point.<\/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\/189\/findFP3_03.png\"> <h3>First Iteration<a name=\"9\"><\/a><\/h3>\r\n   <p>As described in the previous post, I now set the current value of <tt>g2<\/tt> to the new <tt>x<\/tt> value, sliding onto the straight line, and then calculate the next value of <tt>g2<\/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) = g2(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\/189\/findFP3_04.png\"> <h3>Second Iteration<a name=\"10\"><\/a><\/h3>\r\n   <p>Here's the second iteration.<\/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) = g2(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\/189\/findFP3_05.png\"> <h3>Third Iteration<a name=\"11\"><\/a><\/h3>\r\n   <p>And the third.<\/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) = g2(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\/189\/findFP3_06.png\"> <h3>A Few More Rounds<a name=\"12\"><\/a><\/h3>\r\n   <p>Let's do a few more iterations and preallocate the arrays instead of growing them.  And yes, I know in this post I kept overplotting\r\n      lines successively. I didn't want us to get distracted by managing the plots.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">x(9:22) = NaN;\r\ny(9:22) = NaN;\r\n<span style=\"color: #0000FF\">for<\/span> n = 9:2:21\r\n    x(n) = y(n-1);\r\n    y(n) = y(n-1);\r\n    x(n+1) = x(n);\r\n    y(n+1) = g2(x(n+1));\r\n<span style=\"color: #0000FF\">end<\/span>\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\/189\/findFP3_07.png\"> <h3>Fixed Point Found!<a name=\"13\"><\/a><\/h3>\r\n   <p>We appear to be converging on a fixed point, since after iterating, the final values for <tt>x<\/tt> and <tt>g2(x)<\/tt> (which is <tt>y<\/tt>) are getting quite close to each other.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">final = [ x(end) y(end) fzsolution]<\/pre><pre style=\"font-style:oblique\">final =\r\n      0.68039      0.68356      0.68233\r\n<\/pre><h3>Wrapping Up the Series of Posts<a name=\"14\"><\/a><\/h3>\r\n   <p>This post completes this series of posts on fixed point iteration.  There is, of course, more that you could do in a class.\r\n       For example, you might explore what characteristics of the functions <tt>g1<\/tt> and <tt>g2<\/tt> made the fixed point iteration strategy behave differently.  Perhaps look at first derivatives?\r\n   <\/p>\r\n   <p>Here's the recap of the series<\/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<\/a><\/li>\r\n         <li><a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=189\">Rooting Around in MATLAB - Part 3 (this post)<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>Is this sort of tutorial relevant to your work, especially if you are teaching?  How about the incorporation of the visualization\r\n      during the exploration?  Let me know your thoughts <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=189#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_37226cf422ce4413a6469228b4360280() {\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='37226cf422ce4413a6469228b4360280 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 37226cf422ce4413a6469228b4360280';\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_37226cf422ce4413a6469228b4360280()\"><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\n37226cf422ce4413a6469228b4360280 ##### SOURCE BEGIN #####\r\n%% Rooting Around in MATLAB - Part 3\r\n% I recently wrote a pair of posts \r\n% (<https:\/\/blogs.mathworks.com\/loren\/?p=187 1> and \r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=188 2>) about finding roots of \r\n% equations, and how you might use MATLAB to explore this topic with\r\n% emphasis on the method 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 Reminder\r\n% Let me restate the example function. Let's start with a simple cubic\r\n% polynomial |f(x)|\r\n%\r\n% $f(x) = x^3 + x - 1$ \r\n%\r\n% which we define using an anonymous function.\r\nf = @(x) x.^3 + x - 1;\r\nfplot(f,[-2 2])\r\ntitle f\r\ngrid on\r\n%% First Fixed Point Iteration Attempt\r\n% Previously we rewrote |f(x)=0| so that |x| was alone on one side of the\r\n% equation.\r\n%\r\n% $x = 1 - x^3$\r\n%\r\ng1 = @(x) 1 - x.^3;\r\n%%\r\n% This function, |g1(x)| did not help find the root between 0 and 1 - every\r\n% step took us further away from the solutions we found with |roots| and\r\n% |fzero|.\r\nfzsolution = fzero(f,0.5)\r\n%% Second Way to Rewrite Equation\r\n% There's another way we can write the equation for the fixed point.\r\n% Remember, we want to rearrange |f(x)=0| into something like |g2(x)=x|.  \r\n% We have already tried |g1 = 1 - x^3|.  We can also rewrite the equation as\r\n%\r\n% $x = 1\/ (x^2 + 1)$\r\ng2 = @(x) 1.\/(x.^2+1);\r\nfplot(g2,[0 1]);\r\nhold on\r\nstraightLine = @(x) x;\r\nfplot(straightLine, [0 1], 'g')\r\nlegend('g2','x','Location','SouthEast')\r\ngrid on\r\naxis equal, axis([0 1 0 1])\r\n%% Let's Iterate\r\n% Following the same idea from the last post, we now iterate, finding\r\n% successive guesses of |x|, computing the value of |g2(x)|, setting this\r\n% value to be the next estimate of |x|, etc.\r\n%% Let's Try to \"Zero\" in on the Solution\r\n% First, we'll select a starting value |x = 0.5|, and |y|, from the\r\n% straight line also at 0.5.\r\nx(1) = 0.5;\r\ny(1) = x(1);\r\nx(2) = x(1);\r\ny(2) = g2(x(2));\r\n%%\r\n% And let's plot the trajectory of the solutions, starting with the first\r\n% point.\r\nplot(x,y,'r',x(1),y(1),'r*')\r\n%% First Iteration\r\n% As described in the previous post, I now set the current value of |g2| to\r\n% the new |x| value, sliding onto the straight line, and then calculate the\r\n% next value of |g2| 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) = g2(x(4));\r\nplot(x,y,'r')\r\n%% Second Iteration\r\n% Here's the second iteration.\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) = g2(x(n+1));\r\nplot(x,y,'r')\r\n%% Third Iteration\r\n% And the third.\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) = g2(x(n+1));\r\nplot(x,y,'r')\r\n%% A Few More Rounds\r\n% Let's do a few more iterations and preallocate the arrays instead of\r\n% growing them.  And yes, I know in this post I kept overplotting lines\r\n% successively. I didn't want us to get distracted by managing the plots.\r\nx(9:22) = NaN;\r\ny(9:22) = NaN;\r\nfor n = 9:2:21\r\n    x(n) = y(n-1);\r\n    y(n) = y(n-1);\r\n    x(n+1) = x(n);\r\n    y(n+1) = g2(x(n+1));\r\nend\r\nplot(x,y,'r')\r\nhold off\r\n%% Fixed Point Found!\r\n% We appear to be converging on a fixed point, since after iterating, the\r\n% final values for |x| and |g2(x)| (which is |y|) are getting quite close\r\n% to each other.\r\nfinal = [ x(end) y(end) fzsolution]\r\n%% Wrapping Up the Series of Posts\r\n% This post completes this series of posts on fixed point iteration.  There\r\n% is, of course, more that you could do in a class.  For example, you might\r\n% explore what characteristics of the functions |g1| and |g2| made the\r\n% fixed point iteration strategy behave differently.  Perhaps look at first\r\n% derivatives?\r\n%%\r\n% Here's the recap of the series\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>\r\n% * <https:\/\/blogs.mathworks.com\/loren\/?p=189 Rooting Around in MATLAB - Part 3 (this post)>\r\n%\r\n% Is this sort of tutorial relevant to your work, especially if you are\r\n% teaching?  How about the incorporation of the visualization during the\r\n% exploration?  Let me know your thoughts\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=189#respond here>.\r\n\r\n\r\n\r\n \r\n\r\n\r\n\r\n\r\n\r\n##### SOURCE END ##### 37226cf422ce4413a6469228b4360280\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      I recently wrote a pair of posts (1 and 2) about finding roots of equations, and how you might use MATLAB to explore this topic with emphasis on the method of fixed point iteration.\r\n   ... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2009\/06\/23\/rooting-around-in-matlab-part-3\/\">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\/189"}],"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=189"}],"version-history":[{"count":0,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/189\/revisions"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=189"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=189"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=189"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}