{"id":6301,"date":"2015-12-11T09:00:58","date_gmt":"2015-12-11T14:00:58","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=6301"},"modified":"2015-12-10T19:46:49","modified_gmt":"2015-12-11T00:46:49","slug":"polynomial-fit-passing-through-specified-points","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2015\/12\/11\/polynomial-fit-passing-through-specified-points\/","title":{"rendered":"Polynomial fit passing through specified points"},"content":{"rendered":"\r\n\r\n<div class=\"content\"><!--introduction--><p><a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/15007\">Jiro<\/a>'s pick this week is <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/54207-polyfix-x-y-n-xfix-yfix-xder-dydx-\"><tt>polyfix<\/tt><\/a> by <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/3934907\">Are Mjaavatten<\/a>.<\/p><p>Have you ever wanted to fit a polynomial to your data and have the line go through some specified points? What about specifying the slope at a certain point? Let's take a look at some options, including Are's entry.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#e2b57899-d390-47a6-b653-c63595c2906e\">Polynomial fitting<\/a><\/li><li><a href=\"#29663c7c-d9df-412e-9844-01787636ae18\">Constrain to go through certain points<\/a><\/li><li><a href=\"#44533b87-5048-4e0b-9fa9-78e691c85797\">Constrain to have a certain derivative<\/a><\/li><li><a href=\"#a460f527-bf0f-4c66-baba-fb4ccbe6a958\"><tt>polyfix<\/tt><\/a><\/li><li><a href=\"#e5cbd006-b67a-49a1-9135-650a248bc971\">Optimization Toolbox<\/a><\/li><li><a href=\"#731b8b35-a96a-4fcd-ad42-7c7d7b9d2f86\">Comments<\/a><\/li><\/ul><\/div><h4>Polynomial fitting<a name=\"e2b57899-d390-47a6-b653-c63595c2906e\"><\/a><\/h4><p>The function <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyfit.html\"><tt>polyfit<\/tt><\/a> lets you fit a polynomial to your data. We'll first create some data.<\/p><pre class=\"codeinput\">t = linspace(0,2,31)';\r\ny = sin(pi*t)+ 0.1*randn(size(t));\r\nplot(t,y,<span class=\"string\">'o'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_01.png\" alt=\"\"> <p>We'll fit a 3rd order polynomial to the data.<\/p><p>$$ y \\approx x_1t^3 + x_2t^2 + x_3t + x_4$$<\/p><pre class=\"codeinput\">x_polyfit = polyfit(t,y,3)\r\ny1 = polyval(x_polyfit,t);\r\nplot(t,y,<span class=\"string\">'o'<\/span>,t,y1)\r\nlegend(<span class=\"string\">'data'<\/span>,<span class=\"string\">'polyfit'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">x_polyfit =\r\n    2.8843   -8.6941    6.0419   -0.2050\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_02.png\" alt=\"\"> <h4>Constrain to go through certain points<a name=\"29663c7c-d9df-412e-9844-01787636ae18\"><\/a><\/h4><p>What if you want this polynomial to go through certain points. Perhaps, you want the curve to cross (0, 0) and (2, 0). This is where Are's entry comes into play. But first, let me talk about a different method. I found <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/94272-how-do-i-constrain-a-fitted-curve-through-specific-points-like-the-origin-in-matlab\">this question<\/a> on <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/\">MATLAB Answers<\/a>. There are several ways to deal with this, and one of them is to use a function like <a href=\"https:\/\/www.mathworks.com\/help\/optim\/ug\/lsqlin.html\"><tt>lsqlin<\/tt><\/a> from <a href=\"https:\/\/www.mathworks.com\/products\/optimization\/\">Optimization Toolbox<\/a>. <tt>lsqlin<\/tt> solves the following least-squares curve fitting problem.<\/p><p>$$\\min_{x}\\frac{1}{2}||C\\cdot x - d||_2^2 \\quad \\rm {such \\; that} \\quad \\Bigg \\{  \\begin{array}{l} A \\cdot x \\leq b, \\\\ Aeq \\cdot x = beq, \\\\ lb \\leq x \\leq ub \\end{array} $$<\/p><p>For a 3rd order polynomial, $C$ and $d$ are defined as the following.<\/p><p>$$C = \\left[ t^3 \\quad t^2 \\quad t \\quad 1 \\right],\\quad d = y$$<\/p><p>$x$ represents the coefficients for the polynomial terms.<\/p><p>$$x : \\left[ \\begin{array}{c} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\end{array} \\right]$$<\/p><pre class=\"codeinput\">C = [t.^3 t.^2 t ones(size(t))];\r\nd = y;\r\n<\/pre><p>There are no inequality constraints, so we set $A$ and $b$ to empty.<\/p><pre class=\"codeinput\">A = [];     <span class=\"comment\">% No inequality constraint<\/span>\r\nb = [];     <span class=\"comment\">% No inequality constraint<\/span>\r\n<\/pre><p>For the equality constraints, we need to set the $y$ values at $t$ = 0 and $t$ = 2 to zero. This means<\/p><p>$$\\begin{array}{rccccc}  x_1(0)^3 + x_2(0)^2 + x_3(0) + x_4 = 0\r\n&amp; \\Rightarrow &amp;\r\n[0\\quad0\\quad0\\quad1] &amp; \\cdot &amp; x = &amp; 0\r\n\\\\\r\nx_1(2)^3 + x_2(2)^2 + x_3(2) + x_4 = 0\r\n&amp; \\Rightarrow &amp;\r\n\\underbrace{[2^3\\quad2^2\\quad2\\quad1]}_{Aeq} &amp; \\cdot &amp; x = &amp; \\underbrace{0}_{beq} \\end{array}$$<\/p><p>Putting these into $Aeq$ and $beq$, we get the following.<\/p><pre class=\"codeinput\">Aeq = [0    0    0    1 ;   <span class=\"comment\">% t = 0<\/span>\r\n       2^3  2^2  2    1];   <span class=\"comment\">% t = 2<\/span>\r\nbeq = [0 ;                  <span class=\"comment\">% y = 0 (when t = 0)<\/span>\r\n       0];                  <span class=\"comment\">% y = 0 (when t = 2)<\/span>\r\n<\/pre><p>We then call <tt>lsqlin<\/tt> to solve for the polynomial coefficients.<\/p><pre class=\"codeinput\">x_lsqlin1 = lsqlin(C,d,A,b,Aeq,beq)\r\n<\/pre><pre class=\"codeoutput\">Warning: The trust-region-reflective algorithm can handle bound constraints\r\nonly;\r\n    using active-set algorithm instead. \r\nOptimization terminated.\r\nx_lsqlin1 =\r\n    2.5524\r\n   -7.6807\r\n    5.1518\r\n         0\r\n<\/pre><pre class=\"codeinput\">plot(t,y,<span class=\"string\">'o'<\/span>,t,C*x_lsqlin1)\r\nlegend(<span class=\"string\">'data'<\/span>,<span class=\"string\">'lsqlin'<\/span>)\r\ngrid <span class=\"string\">on<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_03.png\" alt=\"\"> <h4>Constrain to have a certain derivative<a name=\"44533b87-5048-4e0b-9fa9-78e691c85797\"><\/a><\/h4><p>Let's add more constraint for the curve to go through (0.5,1) and (1.5,-1). Also, we want the derivatives at those points to be zero. That's additional 4 constraints. With so many constraints, we won't be able to satisfy them with a 3rd order polynomial. Let's try a 7th order.<\/p><p>$$ y \\approx x_1t^7 + x_2t^6 + x_3t^5 + x_4t^4 + x_5t^3 + x_6t^2 + x_7t + x_8$$<\/p><pre class=\"codeinput\">C = [t.^7 t.^6 t.^5 t.^4 t.^3 t.^2 t ones(size(t))];\r\n<\/pre><p>The derivative for the 7th order polynomial is defined as<\/p><p>$$\\frac{dy}{dt} \\approx 7x_1t^6 + 6x_2t^5 + 5x_3t^4 + 4x_4t^3 + 3x_5t^2 + 2x_6t + x_7$$<\/p><p>For example, the constraint <em>\"the derivative at $t$ = 0.5 is zero\"<\/em> can be represented as<\/p><p>$$ \\begin{array}{c} { 7x_1(0.5)^6 + 6x_2(0.5)^5 + 5x_3(0.5)^4 + 4x_4(0.5)^3 + 3x_5(0.5)^2 + 2x_6(0.5) + x_7 = 0 }\r\n\\\\ \\\\ { \\Downarrow } \\\\ \\\\\r\n{ \\underbrace{\\left[ 7(0.5)^6\\quad6(0.5)^5\\quad5(0.5)^4\\quad4(0.5)^3\\quad3(0.5)^2\\quad2(0.5)\\quad1\\quad0\\right]}_{Aeq}\\cdot x = \\underbrace{0}_{beq} } \\end{array}$$<\/p><pre class=\"codeinput\">Aeq = [0       0       0       0       0       0     0     1 ;   <span class=\"comment\">% t = 0<\/span>\r\n       2^7     2^6     2^5     2^4     2^3     2^2   2     1 ;   <span class=\"comment\">% t = 2<\/span>\r\n       0.5^7   0.5^6   0.5^5   0.5^4   0.5^3   0.5^2 0.5   1 ;   <span class=\"comment\">% t = 0.5<\/span>\r\n       1.5^7   1.5^6   1.5^5   1.5^4   1.5^3   1.5^2 1.5   1 ;   <span class=\"comment\">% t = 1.5<\/span>\r\n       7*0.5^6 6*0.5^5 5*0.5^4 4*0.5^3 3*0.5^2 2*0.5 1     0 ;   <span class=\"comment\">% t = 0.5 (for dy\/dt)<\/span>\r\n       7*1.5^6 6*1.5^5 5*1.5^4 4*1.5^3 3*1.5^2 2*1.5 1     0];   <span class=\"comment\">% t = 1.5 (for dy\/dt)<\/span>\r\nbeq = [ 0 ;     <span class=\"comment\">% y = 0 (when t = 0)<\/span>\r\n        0 ;     <span class=\"comment\">% y = 0 (when t = 2)<\/span>\r\n        1 ;     <span class=\"comment\">% y = 1 (when t = 0.5)<\/span>\r\n       -1 ;     <span class=\"comment\">% y = -1 (when t = 1.5)<\/span>\r\n        0 ;     <span class=\"comment\">% dy\/dt = 0 (when t = 0.5)<\/span>\r\n        0];     <span class=\"comment\">% dy\/dt = 0 (when t = 1.5)<\/span>\r\n<\/pre><p>Let's fit and see if we've accomplished our goal.<\/p><pre class=\"codeinput\">x_lsqlin2 = lsqlin(C,d,A,b,Aeq,beq)\r\nplot(t,y,<span class=\"string\">'o'<\/span>,t,C*x_lsqlin2)\r\nlegend(<span class=\"string\">'data'<\/span>,<span class=\"string\">'lsqlin (with derivative constraint)'<\/span>)\r\ngrid <span class=\"string\">on<\/span>\r\nset(gca,<span class=\"string\">'XTick'<\/span>,0:0.5:2)   <span class=\"comment\">% Adjust tick marks to show points of interest<\/span>\r\n<\/pre><pre class=\"codeoutput\">Warning: The trust-region-reflective algorithm can handle bound constraints\r\nonly;\r\n    using active-set algorithm instead. \r\nOptimization terminated.\r\nx_lsqlin2 =\r\n   -0.2925\r\n    2.1030\r\n   -7.8142\r\n   17.6810\r\n  -19.6794\r\n    5.7234\r\n    2.2753\r\n    0.0000\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_04.png\" alt=\"\"> <h4><tt>polyfix<\/tt><a name=\"a460f527-bf0f-4c66-baba-fb4ccbe6a958\"><\/a><\/h4><p>In the MATLAB Answers <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/answers\/94272-how-do-i-constrain-a-fitted-curve-through-specific-points-like-the-origin-in-matlab\">post<\/a> I mentioned above, Are actually posted a response mentioning <tt>polyfix<\/tt>. This entry achieves the goal of performing a polynomial fit with constraints to pass through specific points with specific derivatives. Let's solve the same problem using <tt>polyfix<\/tt>, 7th order polynomial to fit through (0,0), (2,0), (0.5,1), (1.5,-1) and derivative of zero at $t$ = 0.5 and $t$ = 1.5.<\/p><pre class=\"codeinput\">x_polyfix = polyfix(t,y,7,[0 2 0.5 1.5],[0 0 1 -1],[0.5 1.5],[0 0])\r\nplot(t,y,<span class=\"string\">'o'<\/span>,t,polyval(x_polyfix,t))\r\nlegend(<span class=\"string\">'data'<\/span>,<span class=\"string\">'polyfix'<\/span>)\r\ngrid <span class=\"string\">on<\/span>\r\nset(gca,<span class=\"string\">'XTick'<\/span>,0:0.5:2)   <span class=\"comment\">% Adjust tick marks to show points of interest<\/span>\r\n<\/pre><pre class=\"codeoutput\">x_polyfix =\r\n  Columns 1 through 7\r\n   -0.2925    2.1030   -7.8142   17.6810  -19.6794    5.7234    2.2753\r\n  Column 8\r\n         0\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_05.png\" alt=\"\"> <p>If we compare this result with that from <tt>lsqlin<\/tt>, we see that it's essentially identical.<\/p><pre class=\"codeinput\">norm(x_lsqlin2-x_polyfix')\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n   2.1194e-10\r\n<\/pre><h4>Optimization Toolbox<a name=\"e5cbd006-b67a-49a1-9135-650a248bc971\"><\/a><\/h4><p>Are's <tt>polyfit<\/tt> is great for performing polynomial fits with constraints around passing points. If you need other sophisticated constraints, you would want to check out <a href=\"https:\/\/www.mathworks.com\/products\/optimization\/\">Optimization Toolbox<\/a>.<\/p><h4>Comments<a name=\"731b8b35-a96a-4fcd-ad42-7c7d7b9d2f86\"><\/a><\/h4><p>Give this a try and let us know what you think <a href=\"https:\/\/blogs.mathworks.com\/pick\/?p=6301#respond\">here<\/a> or leave a <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/54207-polyfix-x-y-n-xfix-yfix-xder-dydx-#comments\">comment<\/a> for Are.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_fe710d4667914c7085ddc298703cc95c() {\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='fe710d4667914c7085ddc298703cc95c ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' fe710d4667914c7085ddc298703cc95c';\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        copyright = 'Copyright 2015 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 copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\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     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_fe710d4667914c7085ddc298703cc95c()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2015b<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; R2015b<br><\/p><\/div><!--\r\nfe710d4667914c7085ddc298703cc95c ##### SOURCE BEGIN #####\r\n%%\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/15007\r\n% Jiro>'s pick this week is\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/54207-polyfix-x-y-n-xfix-yfix-xder-dydx- |polyfix|> by\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/3934907 Are\r\n% Mjaavatten>.\r\n%\r\n% Have you ever wanted to fit a polynomial to your data and have the line\r\n% go through some specified points? What about specifying the slope at a\r\n% certain point? Let's take a look at some options, including Are's entry.\r\n\r\n%% Polynomial fitting\r\n% The function <https:\/\/www.mathworks.com\/help\/matlab\/ref\/polyfit.html\r\n% |polyfit|> lets you fit a polynomial to your data. We'll first create\r\n% some data.\r\n\r\nt = linspace(0,2,31)';\r\ny = sin(pi*t)+ 0.1*randn(size(t));\r\nplot(t,y,'o')\r\n\r\n%%\r\n% We'll fit a 3rd order polynomial to the data.\r\n%\r\n% $$ y \\approx x_1t^3 + x_2t^2 + x_3t + x_4$$\r\n\r\nx_polyfit = polyfit(t,y,3)\r\ny1 = polyval(x_polyfit,t);\r\nplot(t,y,'o',t,y1)\r\nlegend('data','polyfit')\r\n\r\n%% Constrain to go through certain points\r\n% What if you want this polynomial to go through certain points. Perhaps,\r\n% you want the curve to cross (0, 0) and (2, 0). This is where Are's entry\r\n% comes into play. But first, let me talk about a different method. I found\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/94272-how-do-i-constrain-a-fitted-curve-through-specific-points-like-the-origin-in-matlab this question> on\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/ MATLAB Answers>. There\r\n% are several ways to deal with this, and one of them is to use a function\r\n% like <https:\/\/www.mathworks.com\/help\/optim\/ug\/lsqlin.html |lsqlin|> from\r\n% <https:\/\/www.mathworks.com\/products\/optimization\/ Optimization Toolbox>.\r\n% |lsqlin| solves the following least-squares curve fitting problem.\r\n% \r\n% $$\\min_{x}\\frac{1}{2}||C\\cdot x - d||_2^2 \\quad \\rm {such \\; that} \\quad \\Bigg \\{  \\begin{array}{l} A \\cdot x \\leq b, \\\\ Aeq \\cdot x = beq, \\\\ lb \\leq x \\leq ub \\end{array} $$\r\n%\r\n% For a 3rd order polynomial, $C$ and $d$ are defined as the following.\r\n%\r\n% $$C = \\left[ t^3 \\quad t^2 \\quad t \\quad 1 \\right],\\quad d = y$$\r\n%\r\n% $x$ represents the coefficients for the polynomial terms.\r\n%\r\n% $$x : \\left[ \\begin{array}{c} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\end{array} \\right]$$\r\n\r\nC = [t.^3 t.^2 t ones(size(t))];\r\nd = y;\r\n\r\n%%\r\n% There are no inequality constraints, so we set $A$ and $b$ to empty.\r\n\r\nA = [];     % No inequality constraint\r\nb = [];     % No inequality constraint\r\n\r\n%%\r\n% For the equality constraints, we need to set the $y$ values at $t$ = 0\r\n% and $t$ = 2 to zero. This means\r\n%\r\n% $$\\begin{array}{rccccc}  x_1(0)^3 + x_2(0)^2 + x_3(0) + x_4 = 0\r\n% & \\Rightarrow &\r\n% [0\\quad0\\quad0\\quad1] & \\cdot & x = & 0 \r\n% \\\\\r\n% x_1(2)^3 + x_2(2)^2 + x_3(2) + x_4 = 0\r\n% & \\Rightarrow &\r\n% \\underbrace{[2^3\\quad2^2\\quad2\\quad1]}_{Aeq} & \\cdot & x = & \\underbrace{0}_{beq} \\end{array}$$\r\n%\r\n% Putting these into $Aeq$ and $beq$, we get the following.\r\n\r\nAeq = [0    0    0    1 ;   % t = 0\r\n       2^3  2^2  2    1];   % t = 2\r\nbeq = [0 ;                  % y = 0 (when t = 0)\r\n       0];                  % y = 0 (when t = 2)\r\n\r\n%%\r\n% We then call |lsqlin| to solve for the polynomial coefficients.\r\nx_lsqlin1 = lsqlin(C,d,A,b,Aeq,beq)\r\n\r\n%%\r\nplot(t,y,'o',t,C*x_lsqlin1)\r\nlegend('data','lsqlin')\r\ngrid on\r\n\r\n%% Constrain to have a certain derivative\r\n% Let's add more constraint for the curve to go through (0.5,1) and\r\n% (1.5,-1). Also, we want the derivatives at those points to be zero.\r\n% That's additional 4 constraints. With so many constraints, we won't be\r\n% able to satisfy them with a 3rd order polynomial. Let's try a 7th order.\r\n%\r\n% $$ y \\approx x_1t^7 + x_2t^6 + x_3t^5 + x_4t^4 + x_5t^3 + x_6t^2 + x_7t + x_8$$\r\n\r\nC = [t.^7 t.^6 t.^5 t.^4 t.^3 t.^2 t ones(size(t))];\r\n\r\n%%\r\n% The derivative for the 7th order polynomial is defined as\r\n%\r\n% $$\\frac{dy}{dt} \\approx 7x_1t^6 + 6x_2t^5 + 5x_3t^4 + 4x_4t^3 + 3x_5t^2 + 2x_6t + x_7$$\r\n%\r\n% For example, the constraint _\"the derivative at $t$ = 0.5 is zero\"_ can be\r\n% represented as\r\n%\r\n% $$ \\begin{array}{c} { 7x_1(0.5)^6 + 6x_2(0.5)^5 + 5x_3(0.5)^4 + 4x_4(0.5)^3 + 3x_5(0.5)^2 + 2x_6(0.5) + x_7 = 0 }\r\n% \\\\ \\\\ { \\Downarrow } \\\\ \\\\\r\n% { \\underbrace{\\left[ 7(0.5)^6\\quad6(0.5)^5\\quad5(0.5)^4\\quad4(0.5)^3\\quad3(0.5)^2\\quad2(0.5)\\quad1\\quad0\\right]}_{Aeq}\\cdot x = \\underbrace{0}_{beq} } \\end{array}$$\r\n%\r\n\r\nAeq = [0       0       0       0       0       0     0     1 ;   % t = 0\r\n       2^7     2^6     2^5     2^4     2^3     2^2   2     1 ;   % t = 2\r\n       0.5^7   0.5^6   0.5^5   0.5^4   0.5^3   0.5^2 0.5   1 ;   % t = 0.5\r\n       1.5^7   1.5^6   1.5^5   1.5^4   1.5^3   1.5^2 1.5   1 ;   % t = 1.5\r\n       7*0.5^6 6*0.5^5 5*0.5^4 4*0.5^3 3*0.5^2 2*0.5 1     0 ;   % t = 0.5 (for dy\/dt)\r\n       7*1.5^6 6*1.5^5 5*1.5^4 4*1.5^3 3*1.5^2 2*1.5 1     0];   % t = 1.5 (for dy\/dt)\r\nbeq = [ 0 ;     % y = 0 (when t = 0)\r\n        0 ;     % y = 0 (when t = 2)\r\n        1 ;     % y = 1 (when t = 0.5)\r\n       -1 ;     % y = -1 (when t = 1.5)\r\n        0 ;     % dy\/dt = 0 (when t = 0.5)\r\n        0];     % dy\/dt = 0 (when t = 1.5)\r\n\r\n%%\r\n% Let's fit and see if we've accomplished our goal.\r\n\r\nx_lsqlin2 = lsqlin(C,d,A,b,Aeq,beq)\r\nplot(t,y,'o',t,C*x_lsqlin2)\r\nlegend('data','lsqlin (with derivative constraint)')\r\ngrid on\r\nset(gca,'XTick',0:0.5:2)   % Adjust tick marks to show points of interest\r\n\r\n%% |polyfix|\r\n% In the MATLAB Answers\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/answers\/94272-how-do-i-constrain-a-fitted-curve-through-specific-points-like-the-origin-in-matlab post> I mentioned\r\n% above, Are actually posted a response mentioning |polyfix|. This entry\r\n% achieves the goal of performing a polynomial fit with constraints to pass\r\n% through specific points with specific derivatives. Let's solve the same\r\n% problem using |polyfix|, 7th order polynomial to fit through (0,0),\r\n% (2,0), (0.5,1), (1.5,-1) and derivative of zero at $t$ = 0.5 and $t$ =\r\n% 1.5.\r\n\r\nx_polyfix = polyfix(t,y,7,[0 2 0.5 1.5],[0 0 1 -1],[0.5 1.5],[0 0])\r\nplot(t,y,'o',t,polyval(x_polyfix,t))\r\nlegend('data','polyfix')\r\ngrid on\r\nset(gca,'XTick',0:0.5:2)   % Adjust tick marks to show points of interest\r\n\r\n%%\r\n% If we compare this result with that from |lsqlin|, we see that it's\r\n% essentially identical.\r\nnorm(x_lsqlin2-x_polyfix')\r\n\r\n%% Optimization Toolbox\r\n% Are's |polyfit| is great for performing polynomial fits with constraints\r\n% around passing points. If you need other sophisticated constraints, you\r\n% would want to check out <https:\/\/www.mathworks.com\/products\/optimization\/\r\n% Optimization Toolbox>.\r\n\r\n%% Comments\r\n%\r\n% Give this a try and let us know what you think\r\n% <https:\/\/blogs.mathworks.com\/pick\/?p=6301#respond here> or leave a\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/54207-polyfix-x-y-n-xfix-yfix-xder-dydx-#comments\r\n% comment> for Are.\r\n\r\n##### SOURCE END ##### fe710d4667914c7085ddc298703cc95c\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/pick\/jiro\/potw_polyfix\/potw_polyfix_01.png\" onError=\"this.style.display ='none';\" \/><\/div><p>\r\n\r\nJiro's pick this week is polyfix by Are Mjaavatten.Have you ever wanted to fit a polynomial to your data and have the line go through some specified points? What about specifying the slope at a... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/pick\/2015\/12\/11\/polynomial-fit-passing-through-specified-points\/\">read more >><\/a><\/p>","protected":false},"author":35,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/6301"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/users\/35"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/comments?post=6301"}],"version-history":[{"count":9,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/6301\/revisions"}],"predecessor-version":[{"id":6310,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/posts\/6301\/revisions\/6310"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/media?parent=6301"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/categories?post=6301"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/pick\/wp-json\/wp\/v2\/tags?post=6301"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}