{"id":4707,"date":"2019-04-29T12:00:11","date_gmt":"2019-04-29T17:00:11","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=4707"},"modified":"2019-04-22T21:23:27","modified_gmt":"2019-04-23T02:23:27","slug":"makima-piecewise-cubic-interpolation","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2019\/04\/29\/makima-piecewise-cubic-interpolation\/","title":{"rendered":"Makima Piecewise Cubic Interpolation"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p><i>This post is by my colleague <b>Cosmin Ionita<\/b>.<\/i><\/p><p>The <tt>'makima'<\/tt> cubic interpolation method was recently introduced in MATLAB&reg; in the R2017b release as a new option in <tt>interp1<\/tt>, <tt>interp2<\/tt>, <tt>interp3<\/tt>, <tt>interpn<\/tt>, and <tt>griddedInterpolant<\/tt>. Its implementation is not user visible; thus, we have been receiving inquiries from our users about the specifics of this new cubic method.<\/p><p>In the following, we address our users' inquiries by answering these questions:<\/p><div><ol><li>What is the <tt>'makima'<\/tt> formula?<\/li><li>How does <tt>'makima'<\/tt> compare with MATLAB's other cubic methods?<\/li><\/ol><\/div><p>In a nutshell, <tt>'makima'<\/tt> is short for <i>modified Akima piecewise cubic Hermite interpolation<\/i>. It represents a MATLAB-specific modification of Akima's derivative formula and has the following key properties:<\/p><div><ul><li>It produces undulations which find a nice middle ground between <tt>'spline'<\/tt> and <tt>'pchip'<\/tt>.<\/li><li>It is a local cubic interpolant which generalizes to 2-D grids and higher-dimensional n-D grids.<\/li><li>It increases the robustness of Akima's formula in the edge case of equal side slopes.<\/li><li>It eliminates a special type of overshoot arising when the data is constant for more than two consecutive nodes.<\/li><\/ul><\/div><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#855d8231-5957-4d70-9cb0-874725181ef5\">Akima piecewise cubic Hermite interpolation<\/a><\/li><li><a href=\"#d9a97978-0b09-4a1f-a6a5-504d088631d0\">Modified Akima interpolation -- <tt>'makima'<\/tt><\/a><\/li><li><a href=\"#3b78d49a-0e90-4c1d-9f16-49a23a175132\">Modified Akima formula<\/a><\/li><li><a href=\"#76483eda-efc1-4ad2-bc83-b6b7884157a7\">Generalization to n-D grids<\/a><\/li><li><a href=\"#03107fff-6ad4-4aae-a5a6-fba43e59a4f0\">Summary<\/a><\/li><li><a href=\"#472581be-c4e7-4ce6-b703-70844871b37b\"><tt>code<\/tt><\/a><\/li><li><a href=\"#729b8cf3-89a3-4bcf-812b-2dcb669e967a\"><tt>compareCubicPlots<\/tt><\/a><\/li><li><a href=\"#53876a13-2ec0-4332-b281-370e04da1d6b\"><tt>akima<\/tt><\/a><\/li><li><a href=\"#df0aa10a-1002-47fd-998d-ef30d95b3554\"><tt>akimaSlopes<\/tt><\/a><\/li><li><a href=\"#f0756a08-9e13-481f-b075-3407e814dee6\"><tt>makima<\/tt><\/a><\/li><li><a href=\"#e727bbac-27fa-423a-bb24-45daf57c1101\"><tt>makimaSlopes<\/tt><\/a><\/li><\/ul><\/div><h4>Akima piecewise cubic Hermite interpolation<a name=\"855d8231-5957-4d70-9cb0-874725181ef5\"><\/a><\/h4><p>For each interval $[x_i~x_{i+1})$ of an input data set of nodes $x$ and values $v$, piecewise cubic Hermite interpolation finds a cubic polynomial which not only interpolates the given data values $v_i$ and $v_{i+1}$ at the interval's nodes $x_i$ and $x_{i+1}$, but also has specific derivatives $d_i$ and $d_{i+1}$ at $x_i$ and $x_{i+1}$. For more details we refer to the <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/interp.pdf\">interpolation chapter<\/a> in Cleve's <i>Numerical Computing with MATLAB<\/i> textbook.<\/p><p>The key to cubic Hermite interpolation is the choice of derivatives $d_i$.<\/p><p>In his 1970 <a href=\"https:\/\/dl.acm.org\/citation.cfm?id=321609\">paper<\/a>, Hiroshi Akima introduced a derivative formula which <i>avoids excessive local undulations:<\/i><\/p><p><i>H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.<\/i><\/p><p>Let $\\delta_i =\\left(v_{i+1} -v_i \\right)\/\\left(x_{i+1} -x_i \\right)$ be the slope of the interval $[x_i~x_{i+1})$. Akima's derivative at $x_i$ is defined as:<\/p><p>$$d_i = \\frac{|\\delta_{i+1} - \\delta_i|\\delta_{i-1} +\r\n              |\\delta_{i-1} - \\delta_{i-2}|\\delta_i}\r\n        {|\\delta_{i+1} - \\delta_i| + |\\delta_{i-1} - \\delta_{i-2}|},$$<\/p><p>and represents a weighted average between the slopes $\\delta_{i-1}$ and $\\delta_i$ of the intervals $[x_{i-1}~x_i)$ and $[x_i~x_{i+1})$:<\/p><p>$$d_i = \\frac{w_1}{w_1+w_2}\\delta_{i-1} + \\frac{w_2}{w_1+w_2}\\delta_i \\,,\r\n  \\quad\\quad w_1 = |\\delta_{i+1} - \\delta_i|,\r\n  \\quad\\quad w_2 = |\\delta_{i-1} - \\delta_{i-2}|.$$<\/p><p>Notice that Akima's derivative at $x_i$ is computed locally from the five points $x_{i-2}$, $x_{i-1}$, $x_i$, $x_{i+1}$, and $x_{i+2}$. For the end points $x_1$ and $x_n$, it requires the slopes $\\delta_{-1}, \\delta_0$ and $\\delta_n ,\\delta_{n+1}$. Since these slopes are not available in the input data, Akima proposed using quadratic extrapolation to compute them as $\\delta_0 = 2\\delta_1 -\\delta_2$, $\\delta_{-1} = 2\\delta_0 -\\delta_1$ and $\\delta_n = 2\\delta_{n-1} -\\delta_{n-2}$, $\\delta_{n+1} = 2\\delta_n -\\delta_{n-1}$.<\/p><p>But wait!<\/p><p>MATLAB already has two cubic Hermite interpolation methods (see Cleve's blog <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/16\/splines-and-pchips\/\">Splines and Pchips<\/a>):<\/p><div><ul><li><tt>'spline'<\/tt> computes the derivatives by imposing the constraint of continuous second derivatives (this guarantees a very smooth interpolation result),<\/li><li><tt>'pchip'<\/tt> computes the derivatives by imposing local monotonicity in each interval $[x_i~x_{i+1})$ (this preserves the shape of the data).<\/li><\/ul><\/div><p>How does Akima's formula compare with <tt>'spline'<\/tt> and <tt>'pchip'<\/tt>?<\/p><p>Akima's formula finds a nice middle ground between <tt>'spline'<\/tt> and <tt>'pchip'<\/tt>:<\/p><div><ul><li>Its undulations (or wiggles) have smaller amplitudes than <tt>'spline'<\/tt>.<\/li><li>It is not as aggressive as <tt>'pchip'<\/tt> in reducing the undulations.<\/li><\/ul><\/div><p>Here's a representative example comparing these three cubic Hermite interpolants:<\/p><pre class=\"codeinput\">    x1 = [1 2 3 4 5 5.5 7 8 9 9.5 10];\r\n    v1 = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];\r\n    xq1 = 0.75:0.05:10.25;\r\n    compareCubicPlots(x1,v1,xq1,false,<span class=\"string\">'NE'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_01.png\" alt=\"\"> <h4>Modified Akima interpolation -- <tt>'makima'<\/tt><a name=\"d9a97978-0b09-4a1f-a6a5-504d088631d0\"><\/a><\/h4><p>Akima's formula certainly produces nice results. However, we are not yet ready to fully adopt it.<\/p><p><b>Edge case of equal side slopes<\/b><\/p><p>When the lower slopes are equal, $\\delta_{i-2} =\\delta_{i-1}$, and the upper slopes are equal, $\\delta_i =\\delta_{i+1}$, both the numerator and denominator become equal to 0 and Akima's formula returns <tt>NaN<\/tt>. Akima himself recognized this problem and proposed taking the average of the lower and upper slopes for this edge case: $d_i =\\left(\\delta_{i-1} +\\delta_i \\right)\/2$.<\/p><p>Let's try this fix for the following data set where the intervals $[3~4)$ and $[4~5)$ have equal slopes $\\delta_3 =\\delta_4 =1$ and the intervals $[5~6)$ and $[6~7)$ also have equal slopes $\\delta_5 =\\delta_6 =0$:<\/p><pre class=\"codeinput\">    x2 = 1:8;\r\n    v2 = [-1 -1 -1 0 1 1 1 1];\r\n    xq2 = 0.75:0.05:8.25;\r\n    compareCubicPlots(x2,v2,xq2,false,<span class=\"string\">'SE'<\/span>)\r\n    ylim([-1.2 1.2])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_02.png\" alt=\"\"> <p>In this case, Akima's <tt>NaN<\/tt> derivative at $x_5=5$ gets replaced with the average slope $d_5 =\\left(\\delta_4 +\\delta_5 \\right)\/2=0.5$.<\/p><p>But we are still not ready!<\/p><p>When should we switch from Akima's formula to the averaging formula? What if the slopes are <i>almost<\/i> equal?<\/p><p>Let's change the above example such that $\\delta_5 =\\varepsilon$ and $\\delta_6 =-\\varepsilon$ by setting $v_6 =1+\\varepsilon$ , where $\\varepsilon =2^{-52}$ represents <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/eps.html\">machine epsilon<\/a>. We can now compare Akima's full formula with the averaging formula at $d_5$ for two data sets that differ only by $\\varepsilon$:<\/p><pre class=\"codeinput\"><span class=\"comment\">% v2eps: new data set which differs only by an eps in v2(6)<\/span>\r\n    v2eps = v2; v2eps(6) = v2eps(6) + eps;\r\n\r\n    plot(x2,v2,<span class=\"string\">'ko'<\/span>,<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"string\">'MarkerSize'<\/span>,10,<span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'Input data v2'<\/span>)\r\n    hold <span class=\"string\">on<\/span>\r\n    plot(xq2,akima(x2,v2,xq2),<span class=\"string\">'Color'<\/span>,[1 2\/3 0],<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'v2(6) = 1 uses averaging formula'<\/span>)\r\n    plot(xq2,akima(x2,v2eps,xq2),<span class=\"string\">'-.'<\/span>,<span class=\"string\">'Color'<\/span>,[1 2\/3 0],<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'v2(6) = 1+eps uses Akima''s formula'<\/span>)\r\n    hold <span class=\"string\">off<\/span>\r\n    legend(<span class=\"string\">'Location'<\/span>,<span class=\"string\">'SE'<\/span>)\r\n    title(<span class=\"string\">'Akima interpolant for (almost) equal side slopes'<\/span> )\r\n    ylim([-1.2 1.2])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_03.png\" alt=\"\"> <p>The averaging formula (solid line) returns the derivative $d_5=0.5$ while Akima's formula (dashed line) returns $d_5=1$. Therefore, there is an unwanted non-negligible difference at around $x_5=5$ between the two Akima interpolants, even though the two underlying data sets differ only by $\\varepsilon$ at $v_6$.<\/p><p><b>Requirement (1)<\/b><\/p><p>We do not want our Akima interpolant to switch to a different formula for edge cases.<\/p><p><b>Overshoot if data is constant for more than two consecutive nodes<\/b><\/p><p>The previous example was strategically chosen to also reveal another property of the Akima interpolant: if the data is constant for more than two consecutive nodes, like $v_5=v_6=v_7=1$ in the previous example, then the Akima interpolant may produce an overshoot, namely, the interpolant's undulation in the $[5~6)$ interval above.<\/p><p>This special type of overshoot is not desirable in many engineering applications.<\/p><p>The same concerns hold for the undershoot present in the $[2~3)$ interval above.<\/p><p><b>Requirement (2)<\/b><\/p><p>We do not want our Akima interpolant to produce overshoot (or undershoot) when the data is constant for more than two consecutive nodes.<\/p><h4>Modified Akima formula<a name=\"3b78d49a-0e90-4c1d-9f16-49a23a175132\"><\/a><\/h4><p>In MATLAB, <tt>'makima'<\/tt> cubic Hermite interpolation addresses requirements (1) and (2) outlined above.<\/p><p>To eliminate overshoot and avoid edge cases of both numerator and denominator being equal to 0, we modify Akima's derivative formula by tweaking the weights $w_1$ and $w_2$ of the slopes $\\delta_{i-1}$ and $\\delta_i$:<\/p><p>$$\r\nd_i = \\frac{w_1}{w_1+w_2}\\delta_{i-1} + \\frac{w_2}{w_1+w_2}\\delta_i \\,,\r\n\\quad\\quad w_1 = |\\delta_{i+1} - \\delta_i| + |\\delta_{i+1} + \\delta_i|\/2,\r\n\\quad\\quad\r\nw_2 = |\\delta_{i-1} - \\delta_{i-2}| + |\\delta_{i-1} + \\delta_{i-2}| \/ 2.\r\n$$<\/p><p><i><b>This modified formula represents the <tt>'makima'<\/tt> derivative formula used in MATLAB:<\/b><\/i><\/p><div><ul><li>The addition of the $\\left|\\delta_{i+1} +\\delta_i \\right|\/2$ and $\\left|\\delta_{i-1} +\\delta_{i-2} \\right|\/2$ terms forces $d_i =0$ when $\\delta_i =\\delta_{i+1} =0$, i.e., $d_i =0$ when $v_i =v_{i+1} =v_{i+2}$. Therefore, it eliminates overshoot when the data is constant for more than two consecutive nodes.<\/li><li>These new terms also address the edge case of equal side slopes discussed above, $\\delta_{i-2} =\\delta_{i-1}$ and $\\delta_i =\\delta_{i+1}$. However, there is one case which slips through: if the data is constant $v_{i-2} = v_{i-1} = v_i = v_{i+1} = v_{i+2}$, then the four slopes are all equal to zero and we get $d_i = \\mathrm{NaN}$. For this special case of constant data, we set $d_i =0$.<\/li><\/ul><\/div><p>Let's try the <tt>'makima'<\/tt> formula on the above overshoot example:<\/p><pre class=\"codeinput\">    compareCubicPlots(x2,v2,xq2,true,<span class=\"string\">'SE'<\/span>)\r\n    ylim([-1.2 1.2])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_04.png\" alt=\"\"> <p>Indeed, <tt>'makima'<\/tt> does not produce an overshoot if the data is constant for more than two nodes ($v_5=v_6=v_7=1$ above).<\/p><p>But what does this mean for the undulations we saw in our first example?<\/p><pre class=\"codeinput\">    compareCubicPlots(x1,v1,xq1,true,<span class=\"string\">'NE'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_05.png\" alt=\"\"> <p>Notice that <tt>'makima'<\/tt> closely follows the result obtained with Akima's formula. In fact, the results are so similar that it is hard to tell them apart on the plot.<\/p><p>Therefore, <tt>'makima'<\/tt> still preserves Akima's desirable properties of being a nice middle ground between <tt>'spline'<\/tt> and <tt>'pchip'<\/tt> in terms of the resulting undulations.<\/p><h4>Generalization to n-D grids<a name=\"76483eda-efc1-4ad2-bc83-b6b7884157a7\"><\/a><\/h4><p>But we are still not done!<\/p><p>Akima's formula and our modified <tt>'makima'<\/tt> formula have another desirable property: they generalize to higher dimensional n-D gridded data. Akima noticed this property in his 1974 <a href=\"https:\/\/dl.acm.org\/citation.cfm?id=360779\">paper<\/a>.<\/p><p><i>H. Akima, \"A Method of Bivariate Interpolation and Smooth Surface Fitting Based on Local Procedures\", Communications of the ACM, v. 17-1, p.18-20, 1974.<\/i><\/p><p>For example, to interpolate on a 2-D grid $\\left(x,y\\right)$, we first apply the <tt>'makima'<\/tt> derivative formula separately along $x$ and $y$ to obtain two directional derivatives for each grid node. Then, we also compute the cross-derivative along $xy$ for each grid node. The cross-derivative formula first computes the 2-D divided differences corresponding to the 2-D grid data and applies the <tt>'makima'<\/tt> derivative formula along $x$ on these differences; then, it takes the result and applies the derivative formula along $y$. The derivatives and cross-derivatives are then plugged in as coefficients of a two-variable cubic Hermite polynomial representing the 2-D interpolant.<\/p><p>The same procedure applies to higher dimensional n-D grids with $n\\ge 2$ and requires computing cross-derivatives along all possible cross-directions. Therefore, <tt>'makima'<\/tt> is supported not only in <tt>interp1<\/tt>, but also in <tt>interp2<\/tt>, <tt>interp3<\/tt>, <tt>interpn<\/tt>, and <tt>griddedInterpolant<\/tt>.<\/p><p>Here's how 2-D <tt>'makima'<\/tt> interpolation compares with 2-D <tt>'spline'<\/tt> interpolation on gridded data generated with the <tt>peaks<\/tt> function:<\/p><pre class=\"codeinput\">    [X1,Y1,V1] = peaks(5);\r\n    [Xq1,Yq1] = meshgrid(-3:0.1:3,-3:0.1:3);\r\n\r\n    Vqs1 = interp2(X1,Y1,V1,Xq1,Yq1,<span class=\"string\">'spline'<\/span>);\r\n    surf(Xq1,Yq1,Vqs1)\r\n    axis <span class=\"string\">tight<\/span>\r\n    title(<span class=\"string\">'2-D ''spline'''<\/span>)\r\n    snapnow\r\n\r\n    Vqm1 = interp2(X1,Y1,V1,Xq1,Yq1,<span class=\"string\">'makima'<\/span>);\r\n    surf(Xq1,Yq1,Vqm1)\r\n    axis <span class=\"string\">tight<\/span>\r\n    title(<span class=\"string\">'2-D ''makima'''<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_06.png\" alt=\"\"> <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_07.png\" alt=\"\"> <p>Notice the smaller undulations (or wiggles) generated by <tt>'makima'<\/tt>.<\/p><p>Finally, let's try <tt>'makima'<\/tt> on an example with a few 2-D peaks where the data has sharp edges and steps:<\/p><pre class=\"codeinput\">    V2 = zeros(10);\r\n    V2(2:5,2:5) = 3\/7; <span class=\"comment\">% First step<\/span>\r\n    V2(6:7,6:7) = [1\/4 1\/5; 1\/5 1\/4]; <span class=\"comment\">% Middle step<\/span>\r\n    V2(8:9,8:9) = 1\/2; <span class=\"comment\">% Last step<\/span>\r\n    V2 = flip(V2,2);\r\n    [Xq2,Yq2] = meshgrid(1:0.2:10,1:0.2:10);\r\n\r\n    Vqs2 = interp2(V2,Xq2,Yq2,<span class=\"string\">'spline'<\/span>);\r\n    surf(Xq2,Yq2,Vqs2)\r\n    axis <span class=\"string\">tight<\/span>\r\n    title(<span class=\"string\">'2-D ''spline'''<\/span>)\r\n    snapnow\r\n\r\n    Vqm2 = interp2(V2,Xq2,Yq2,<span class=\"string\">'makima'<\/span>);\r\n    surf(Xq2,Yq2,Vqm2)\r\n    axis <span class=\"string\">tight<\/span>\r\n    title(<span class=\"string\">'2-D ''makima'''<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_08.png\" alt=\"\"> <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/makima_blog_09.png\" alt=\"\"> <p>The <tt>'makima'<\/tt> result has smaller undulations (or wiggles) than <tt>'spline'.<\/tt><\/p><h4>Summary<a name=\"03107fff-6ad4-4aae-a5a6-fba43e59a4f0\"><\/a><\/h4><p>We conclude with a summary of properties of interest for the three cubic Hermite interpolation methods supported in MATLAB:<\/p><pre class=\"codeinput\">    summary = table( <span class=\"keyword\">...<\/span>\r\n        [<span class=\"string\">\"C2\"<\/span>; <span class=\"string\">\"Uses all points\"<\/span>; <span class=\"string\">\"Not-a-knot condition\"<\/span>; <span class=\"string\">\"Yes\"<\/span>], <span class=\"keyword\">...<\/span>\r\n        [<span class=\"string\">\"C1\"<\/span>; <span class=\"string\">\"Uses 4 points\"<\/span>; <span class=\"string\">\"One-sided slope\"<\/span>; <span class=\"string\">\"No\"<\/span>], <span class=\"keyword\">...<\/span>\r\n        [<span class=\"string\">\"C1\"<\/span>; <span class=\"string\">\"Uses 5 points\"<\/span>; <span class=\"string\">\"Quadratic extrapolation\"<\/span>; <span class=\"string\">\"Yes\"<\/span>], <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'VariableNames'<\/span>, <span class=\"keyword\">...<\/span>\r\n        [<span class=\"string\">\"spline\"<\/span> <span class=\"string\">\"pchip\"<\/span> <span class=\"string\">\"makima\"<\/span>], <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'RowNames'<\/span>, <span class=\"keyword\">...<\/span>\r\n        [<span class=\"string\">\"Continuity\"<\/span>; <span class=\"string\">\"Locality\"<\/span>; <span class=\"string\">\"End points\"<\/span>; <span class=\"string\">\"Supports n-D\"<\/span>]);\r\n    disp(summary)\r\n<\/pre><pre class=\"codeoutput\">                            spline                  pchip                   makima          \r\n                    ______________________    _________________    _________________________\r\n    Continuity      \"C2\"                      \"C1\"                 \"C1\"                     \r\n    Locality        \"Uses all points\"         \"Uses 4 points\"      \"Uses 5 points\"          \r\n    End points      \"Not-a-knot condition\"    \"One-sided slope\"    \"Quadratic extrapolation\"\r\n    Supports n-D    \"Yes\"                     \"No\"                 \"Yes\"                    \r\n<\/pre><h4><tt>code<\/tt><a name=\"472581be-c4e7-4ce6-b703-70844871b37b\"><\/a><\/h4><h4><tt>compareCubicPlots<\/tt><a name=\"729b8cf3-89a3-4bcf-812b-2dcb669e967a\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> compareCubicPlots(x,v,xq,showMakima,legendLocation)\r\n<span class=\"comment\">% Plot 'pchip', 'spline', Akima, and 'makima' interpolation results<\/span>\r\n    vqp = pchip(x,v,xq);      <span class=\"comment\">% same as vq = interp1(x,v,xq,'pchip')<\/span>\r\n    vqs = spline(x,v,xq);     <span class=\"comment\">% same as vq = interp1(x,v,xq,'spline')<\/span>\r\n    vqa = akima(x,v,xq);\r\n    <span class=\"keyword\">if<\/span> showMakima\r\n        vqm = makima(x,v,xq); <span class=\"comment\">% same as vq = interp1(x,v,xq,'makima')<\/span>\r\n    <span class=\"keyword\">end<\/span>\r\n\r\n    plot(x,v,<span class=\"string\">'ko'<\/span>,<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"string\">'MarkerSize'<\/span>,10,<span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'Input data'<\/span>)\r\n    hold <span class=\"string\">on<\/span>\r\n    plot(xq,vqp,<span class=\"string\">'LineWidth'<\/span>,4,<span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'''pchip'''<\/span>)\r\n    plot(xq,vqs,<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'''spline'''<\/span>)\r\n    plot(xq,vqa,<span class=\"string\">'LineWidth'<\/span>,2,<span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'Akima''s formula'<\/span>)\r\n    <span class=\"keyword\">if<\/span> showMakima\r\n        plot(xq,vqm,<span class=\"string\">'--'<\/span>,<span class=\"string\">'Color'<\/span>,[0 3\/4 0],<span class=\"string\">'LineWidth'<\/span>,2, <span class=\"keyword\">...<\/span>\r\n            <span class=\"string\">'DisplayName'<\/span>,<span class=\"string\">'''makima'''<\/span>)\r\n    <span class=\"keyword\">end<\/span>\r\n    hold <span class=\"string\">off<\/span>\r\n    xticks(x(1)-1:x(end)+1)\r\n    legend(<span class=\"string\">'Location'<\/span>,legendLocation)\r\n    title(<span class=\"string\">'Cubic Hermite interpolation'<\/span>)\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4><tt>akima<\/tt><a name=\"53876a13-2ec0-4332-b281-370e04da1d6b\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> vq = akima(x,v,xq)\r\n<span class=\"comment\">%AKIMA   Akima piecewise cubic Hermite interpolation.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   vq = AKIMA(x,v,xq) interpolates the data (x,v) at the query points xq<\/span>\r\n<span class=\"comment\">%   using Akima's cubic Hermite interpolation formula.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   References:<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting<\/span>\r\n<span class=\"comment\">%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   AKIMA vs. PCHIP vs. SPLINE:<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%     - Akima's cubic formula is a middle ground between SPLINE and PCHIP:<\/span>\r\n<span class=\"comment\">%       It has lower-amplitude wiggles than SPLINE, but is not as agressive<\/span>\r\n<span class=\"comment\">%       at reducing the wiggles as PCHIP.<\/span>\r\n<span class=\"comment\">%     - Akima's cubic formula and SPLINE generalize to n-D grids.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   Example: Compare AKIMA with MAKIMA, SPLINE, and PCHIP<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       x = [1 2 3 4 5 5.5 7 8 9 9.5 10];<\/span>\r\n<span class=\"comment\">%       v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];<\/span>\r\n<span class=\"comment\">%       xq = 0.75:0.05:10.25;<\/span>\r\n<span class=\"comment\">%       vqs = spline(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqp = pchip(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqa = akima(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqm = makima(x,v,xq);<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       figure<\/span>\r\n<span class=\"comment\">%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on<\/span>\r\n<span class=\"comment\">%       plot(xq,vqp,'LineWidth',4)<\/span>\r\n<span class=\"comment\">%       plot(xq,vqs,xq,vqa,xq,vqm,'--','LineWidth',2)<\/span>\r\n<span class=\"comment\">%       legend('Data','''pchip''','''spline''','Akima''s formula','''makima''')<\/span>\r\n<span class=\"comment\">%       title('Cubic Hermite interpolation in MATLAB')<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   See also MAKIMA, SPLINE, PCHIP.<\/span>\r\n\r\n\r\n<span class=\"comment\">%   mailto: cosmin.ionita@mathworks.com<\/span>\r\n\r\n    assert(isvector(x) &amp;&amp; isvector(v) &amp;&amp; (numel(x) == numel(v)))\r\n    assert(numel(x) &gt; 2)\r\n    x = x(:).'; v = v(:).'; <span class=\"comment\">% Same shapes as in pchip.m and spline.m<\/span>\r\n\r\n    <span class=\"comment\">% Compute Akima slopes<\/span>\r\n    h = diff(x);\r\n    delta = diff(v).\/h;\r\n    slopes = akimaSlopes(delta);\r\n\r\n    <span class=\"comment\">% Evaluate the piecewise cubic Hermite interpolant<\/span>\r\n    vq = ppval(pwch(x,v,slopes,h,delta),xq);\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4><tt>akimaSlopes<\/tt><a name=\"df0aa10a-1002-47fd-998d-ef30d95b3554\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> s = akimaSlopes(delta)\r\n<span class=\"comment\">% Derivative values for Akima cubic Hermite interpolation<\/span>\r\n\r\n<span class=\"comment\">% Akima's derivative estimate at grid node x(i) requires the four finite<\/span>\r\n<span class=\"comment\">% differences corresponding to the five grid nodes x(i-2:i+2).<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% For boundary grid nodes x(1:2) and x(n-1:n), append finite differences<\/span>\r\n<span class=\"comment\">% which would correspond to x(-1:0) and x(n+1:n+2) by using the following<\/span>\r\n<span class=\"comment\">% uncentered difference formula correspondin to quadratic extrapolation<\/span>\r\n<span class=\"comment\">% using the quadratic polynomial defined by data at x(1:3)<\/span>\r\n<span class=\"comment\">% (section 2.3 in Akima's paper):<\/span>\r\n    n = numel(delta) + 1; <span class=\"comment\">% number of grid nodes x<\/span>\r\n    delta_0  = 2*delta(1)   - delta(2);\r\n    delta_m1 = 2*delta_0    - delta(1);\r\n    delta_n  = 2*delta(n-1) - delta(n-2);\r\n    delta_n1 = 2*delta_n    - delta(n-1);\r\n    delta = [delta_m1 delta_0 delta delta_n delta_n1];\r\n\r\n<span class=\"comment\">% Akima's derivative estimate formula (equation (1) in the paper):<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting<\/span>\r\n<span class=\"comment\">%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i))<\/span>\r\n<span class=\"comment\">%      \/ (|d(i+1)-d(i)|          + |d(i-1)-d(i-2)|)<\/span>\r\n    weights = abs(diff(delta));\r\n    weights1 = weights(1:n);   <span class=\"comment\">% |d(i-1)-d(i-2)|<\/span>\r\n    weights2 = weights(3:end); <span class=\"comment\">% |d(i+1)-d(i)|<\/span>\r\n    delta1 = delta(2:n+1);     <span class=\"comment\">% d(i-1)<\/span>\r\n    delta2 = delta(3:n+2);     <span class=\"comment\">% d(i)<\/span>\r\n\r\n    weights12 = weights1 + weights2;\r\n    s = (weights2.\/weights12) .* delta1 + (weights1.\/weights12) .* delta2;\r\n\r\n<span class=\"comment\">% To avoid 0\/0, Akima proposed to average the divided differences d(i-1)<\/span>\r\n<span class=\"comment\">% and d(i) for the edge case of d(i-2) = d(i-1) and d(i) = d(i+1):<\/span>\r\n    ind = weights1 == 0 &amp; weights2 == 0;\r\n    s(ind) = (delta1(ind) + delta2(ind)) \/ 2;\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4><tt>makima<\/tt><a name=\"f0756a08-9e13-481f-b075-3407e814dee6\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> vq = makima(x,v,xq)\r\n<span class=\"comment\">%MAKIMA   Modified Akima piecewise cubic Hermite interpolation.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       We modify Akima's formula to eliminate overshoot and undershoot<\/span>\r\n<span class=\"comment\">%       when the data is constant for more than two consecutive nodes.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   vq = MAKIMA(x,v,xq) interpolates the data (x,v) at the query points xq<\/span>\r\n<span class=\"comment\">%   using a modified Akima cubic Hermite interpolation formula.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   References:<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting<\/span>\r\n<span class=\"comment\">%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   MAKIMA vs. PCHIP vs. SPLINE:<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%     - MAKIMA is a middle ground between SPLINE and PCHIP:<\/span>\r\n<span class=\"comment\">%       It has lower-amplitude wiggles than SPLINE, but is not as agressive<\/span>\r\n<span class=\"comment\">%       at reducing the wiggles as PCHIP.<\/span>\r\n<span class=\"comment\">%     - MAKIMA and SPLINE generalize to n-D grids.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   Example: No overshoot and undershoot with MAKIMA when the data is<\/span>\r\n<span class=\"comment\">%            constant for more than two consecutive nodes<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       x = 1:7;<\/span>\r\n<span class=\"comment\">%       v = [-1 -1 -1 0 1 1 1];<\/span>\r\n<span class=\"comment\">%       xq = 0.75:0.05:7.25;<\/span>\r\n<span class=\"comment\">%       vqs = spline(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqp = pchip(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqa = akima(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqm = makima(x,v,xq);<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       figure<\/span>\r\n<span class=\"comment\">%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on<\/span>\r\n<span class=\"comment\">%       plot(xq,vqp,'LineWidth',4)<\/span>\r\n<span class=\"comment\">%       plot(xq,vqs,xq,vqa,xq,vqm,'LineWidth',2)<\/span>\r\n<span class=\"comment\">%       legend('Data','''pchip''','''spline''','Akima''s formula',...<\/span>\r\n<span class=\"comment\">%           '''makima''','Location','SE')<\/span>\r\n<span class=\"comment\">%       title('''makima'' has no overshoot in x(1:3) and x(5:7)')<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   Example: Compare MAKIMA with AKIMA, SPLINE, and PCHIP<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       x = [1 2 3 4 5 5.5 7 8 9 9.5 10];<\/span>\r\n<span class=\"comment\">%       v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];<\/span>\r\n<span class=\"comment\">%       xq = 0.75:0.05:10.25;<\/span>\r\n<span class=\"comment\">%       vqs = spline(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqp = pchip(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqa = akima(x,v,xq);<\/span>\r\n<span class=\"comment\">%       vqm = makima(x,v,xq);<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       figure<\/span>\r\n<span class=\"comment\">%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on<\/span>\r\n<span class=\"comment\">%       plot(xq,vqp,'LineWidth',4)<\/span>\r\n<span class=\"comment\">%       plot(xq,vqs,xq,vqa,xq,vqm,'--','LineWidth',2)<\/span>\r\n<span class=\"comment\">%       legend('Data','''pchip''','''spline''','Akima''s formula','''makima''')<\/span>\r\n<span class=\"comment\">%       title('Cubic Hermite interpolation in MATLAB')<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%   See also AKIMA, SPLINE, PCHIP.<\/span>\r\n\r\n\r\n<span class=\"comment\">%   mailto: cosmin.ionita@mathworks.com<\/span>\r\n\r\n    assert(isvector(x) &amp;&amp; isvector(v) &amp;&amp; (numel(x) == numel(v)))\r\n    assert(numel(x) &gt; 2)\r\n    x = x(:).'; v = v(:).'; <span class=\"comment\">% Same shapes as in pchip.m and spline.m<\/span>\r\n\r\n<span class=\"comment\">% Compute modified Akima slopes<\/span>\r\n    h = diff(x);\r\n    delta = diff(v).\/h;\r\n    slopes = makimaSlopes(delta);\r\n\r\n<span class=\"comment\">% Evaluate the piecewise cubic Hermite interpolant<\/span>\r\n    vq = ppval(pwch(x,v,slopes,h,delta),xq);\r\n\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4><tt>makimaSlopes<\/tt><a name=\"e727bbac-27fa-423a-bb24-45daf57c1101\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> s = makimaSlopes(delta)\r\n<span class=\"comment\">% Derivative values for modified Akima cubic Hermite interpolation<\/span>\r\n\r\n<span class=\"comment\">% Akima's derivative estimate at grid node x(i) requires the four finite<\/span>\r\n<span class=\"comment\">% differences corresponding to the five grid nodes x(i-2:i+2).<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% For boundary grid nodes x(1:2) and x(n-1:n), append finite differences<\/span>\r\n<span class=\"comment\">% which would correspond to x(-1:0) and x(n+1:n+2) by using the following<\/span>\r\n<span class=\"comment\">% uncentered difference formula correspondin to quadratic extrapolation<\/span>\r\n<span class=\"comment\">% using the quadratic polynomial defined by data at x(1:3)<\/span>\r\n<span class=\"comment\">% (section 2.3 in Akima's paper):<\/span>\r\n    n = numel(delta) + 1; <span class=\"comment\">% number of grid nodes x<\/span>\r\n    delta_0  = 2*delta(1)   - delta(2);\r\n    delta_m1 = 2*delta_0    - delta(1);\r\n    delta_n  = 2*delta(n-1) - delta(n-2);\r\n    delta_n1 = 2*delta_n    - delta(n-1);\r\n    delta = [delta_m1 delta_0 delta delta_n delta_n1];\r\n\r\n<span class=\"comment\">% Akima's derivative estimate formula (equation (1) in the paper):<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting<\/span>\r\n<span class=\"comment\">%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i))<\/span>\r\n<span class=\"comment\">%      \/ (|d(i+1)-d(i)|          + |d(i-1)-d(i-2)|)<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% To eliminate overshoot and undershoot when the data is constant for more<\/span>\r\n<span class=\"comment\">% than two consecutive nodes, in MATLAB's 'makima' we modify Akima's<\/span>\r\n<span class=\"comment\">% formula by adding an additional averaging term in the weights.<\/span>\r\n<span class=\"comment\">% s(i) = ( (|d(i+1)-d(i)|   + |d(i+1)+d(i)|\/2  ) * d(i-1) +<\/span>\r\n<span class=\"comment\">%          (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|\/2) * d(i)  )<\/span>\r\n<span class=\"comment\">%      \/ ( (|d(i+1)-d(i)|   + |d(i+1)+d(i)|\/2  ) +<\/span>\r\n<span class=\"comment\">%          (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|\/2)<\/span>\r\n    weights = abs(diff(delta)) + abs((delta(1:end-1)+delta(2:end))\/2);\r\n\r\n    weights1 = weights(1:n);   <span class=\"comment\">% |d(i-1)-d(i-2)|<\/span>\r\n    weights2 = weights(3:end); <span class=\"comment\">% |d(i+1)-d(i)|<\/span>\r\n    delta1 = delta(2:n+1);     <span class=\"comment\">% d(i-1)<\/span>\r\n    delta2 = delta(3:n+2);     <span class=\"comment\">% d(i)<\/span>\r\n\r\n    weights12 = weights1 + weights2;\r\n    s = (weights2.\/weights12) .* delta1 + (weights1.\/weights12) .* delta2;\r\n\r\n<span class=\"comment\">% If the data is constant for more than four consecutive nodes, then the<\/span>\r\n<span class=\"comment\">% denominator is zero and the formula produces an unwanted NaN result.<\/span>\r\n<span class=\"comment\">% Replace this NaN with 0.<\/span>\r\n    s(weights12 == 0) = 0;\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><script language=\"JavaScript\"> <!-- \r\n    function grabCode_e5dc3d9c091f41a6acd2477ff3ef9b53() {\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='e5dc3d9c091f41a6acd2477ff3ef9b53 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' e5dc3d9c091f41a6acd2477ff3ef9b53';\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 2019 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\">Copyright 2019 The MathWorks, Inc.<br><a href=\"javascript:grabCode_e5dc3d9c091f41a6acd2477ff3ef9b53()\"><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; R2018b<br><\/p><\/div><!--\r\ne5dc3d9c091f41a6acd2477ff3ef9b53 ##### SOURCE BEGIN #####\r\n%% The Makima Piecewise Cubic Interpolation Method\r\n% _This post is by my colleague *Cosmin Ionita*._\r\n% \r\n% The |'makima'| cubic interpolation method was recently introduced in\r\n% MATLAB(R) in the R2017b release as a new option in |interp1|, |interp2|,\r\n% |interp3|, |interpn|, and |griddedInterpolant|. Its implementation is not\r\n% user visible; thus, we have been receiving inquiries from our users about\r\n% the specifics of this new cubic method.\r\n% \r\n% In the following, we address our users' inquiries by answering these\r\n% questions:\r\n%\r\n% # What is the |'makima'| formula?\r\n% # How does |'makima'| compare with MATLAB's other cubic methods?\r\n%\r\n% In a nutshell, |'makima'| is short for _modified Akima piecewise cubic\r\n% Hermite interpolation_. It represents a MATLAB-specific modification of\r\n% Akima's derivative formula and has the following key properties:\r\n%\r\n% * It produces undulations which find a nice middle ground between\r\n% |'spline'| and |'pchip'|.\r\n% * It is a local cubic interpolant which generalizes to 2-D grids and\r\n% higher-dimensional n-D grids.\r\n% * It increases the robustness of Akima's formula in the edge case of\r\n% equal side slopes.\r\n% * It eliminates a special type of overshoot arising when the data is\r\n% constant for more than two consecutive nodes.\r\n%\r\n%% Akima piecewise cubic Hermite interpolation\r\n% For each interval $[x_i~x_{i+1})$ of an input data set of nodes $x$ and\r\n% values $v$, piecewise cubic Hermite interpolation finds a cubic\r\n% polynomial which not only interpolates the given data values $v_i$ and\r\n% $v_{i+1}$ at the interval's nodes $x_i$ and $x_{i+1}$, but also has\r\n% specific derivatives $d_i$ and $d_{i+1}$ at $x_i$ and $x_{i+1}$. For more\r\n% details we refer to the\r\n% <https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/interp.pdf\r\n% interpolation chapter> in Cleve's _Numerical Computing with MATLAB_\r\n% textbook.\r\n% \r\n% The key to cubic Hermite interpolation is the choice of derivatives\r\n% $d_i$.\r\n% \r\n% In his 1970 <https:\/\/dl.acm.org\/citation.cfm?id=321609 paper>,\r\n% Hiroshi Akima introduced a derivative formula which _avoids \r\n% excessive local undulations:_\r\n%\r\n% _H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting Based \r\n% on Local Procedures\", JACM, v. 17-4, p.589-602, 1970._\r\n%\r\n% Let $\\delta_i =\\left(v_{i+1} -v_i \\right)\/\\left(x_{i+1} -x_i \\right)$ be\r\n% the slope of the interval $[x_i~x_{i+1})$. Akima's derivative at $x_i$ is\r\n% defined as:\r\n% \r\n% $$d_i = \\frac{|\\delta_{i+1} - \\delta_i|\\delta_{i-1} +\r\n%               |\\delta_{i-1} - \\delta_{i-2}|\\delta_i}\r\n%         {|\\delta_{i+1} - \\delta_i| + |\\delta_{i-1} - \\delta_{i-2}|},$$\r\n% \r\n% and represents a weighted average between the slopes $\\delta_{i-1}$ and\r\n% $\\delta_i$ of the intervals $[x_{i-1}~x_i)$ and $[x_i~x_{i+1})$:\r\n% \r\n% $$d_i = \\frac{w_1}{w_1+w_2}\\delta_{i-1} + \\frac{w_2}{w_1+w_2}\\delta_i \\,, \r\n%   \\quad\\quad w_1 = |\\delta_{i+1} - \\delta_i|,\r\n%   \\quad\\quad w_2 = |\\delta_{i-1} - \\delta_{i-2}|.$$\r\n% \r\n% Notice that Akima's derivative at $x_i$ is computed locally from the five\r\n% points $x_{i-2}$, $x_{i-1}$, $x_i$, $x_{i+1}$, and $x_{i+2}$. For the end\r\n% points $x_1$ and $x_n$, it requires the slopes $\\delta_{-1}, \\delta_0$\r\n% and $\\delta_n ,\\delta_{n+1}$. Since these slopes are not available in the\r\n% input data, Akima proposed using quadratic extrapolation to compute them\r\n% as $\\delta_0 = 2\\delta_1 -\\delta_2$, $\\delta_{-1} = 2\\delta_0 -\\delta_1$\r\n% and $\\delta_n = 2\\delta_{n-1} -\\delta_{n-2}$, $\\delta_{n+1} = 2\\delta_n\r\n% -\\delta_{n-1}$.\r\n% \r\n% But wait!\r\n% \r\n% MATLAB already has two cubic Hermite interpolation methods (see Cleve's\r\n% blog <https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/16\/splines-and-pchips\/\r\n% Splines and Pchips>):\r\n%\r\n% * |'spline'| computes the derivatives by imposing the constraint of\r\n% continuous second derivatives (this guarantees a very smooth\r\n% interpolation result),\r\n% * |'pchip'| computes the derivatives by imposing local monotonicity in\r\n% each interval $[x_i~x_{i+1})$ (this preserves the shape of the data).\r\n%\r\n% How does Akima's formula compare with |'spline'| and |'pchip'|?\r\n% \r\n% Akima's formula finds a nice middle ground between |'spline'| and\r\n% |'pchip'|:\r\n%\r\n% * Its undulations (or wiggles) have smaller amplitudes than |'spline'|.\r\n% * It is not as aggressive as |'pchip'| in reducing the undulations.\r\n%\r\n% Here's a representative example comparing these three cubic Hermite\r\n% interpolants:\r\n\r\n    x1 = [1 2 3 4 5 5.5 7 8 9 9.5 10];\r\n    v1 = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];\r\n    xq1 = 0.75:0.05:10.25;\r\n    compareCubicPlots(x1,v1,xq1,false,'NE')\r\n\r\n%% Modified Akima interpolation REPLACE_WITH_DASH_DASH |'makima'|\r\n% Akima's formula certainly produces nice results. However, we are not yet\r\n% ready to fully adopt it.\r\n%\r\n% *Edge case of equal side slopes*\r\n%\r\n% When the lower slopes are equal, $\\delta_{i-2} =\\delta_{i-1}$, and the\r\n% upper slopes are equal, $\\delta_i =\\delta_{i+1}$, both the numerator and\r\n% denominator become equal to 0 and Akima's formula returns |NaN|.\r\n% Akima himself recognized this problem and proposed taking the average of\r\n% the lower and upper slopes for this edge case: $d_i =\\left(\\delta_{i-1}\r\n% +\\delta_i \\right)\/2$.\r\n% \r\n% Let's try this fix for the following data set where the intervals $[3~4)$\r\n% and $[4~5)$ have equal slopes $\\delta_3 =\\delta_4 =1$ and the intervals\r\n% $[5~6)$ and $[6~7)$ also have equal slopes $\\delta_5 =\\delta_6 =0$:\r\n\r\n    x2 = 1:8;\r\n    v2 = [-1 -1 -1 0 1 1 1 1];\r\n    xq2 = 0.75:0.05:8.25;\r\n    compareCubicPlots(x2,v2,xq2,false,'SE')\r\n    ylim([-1.2 1.2])\r\n    \r\n%%\r\n% In this case, Akima's |NaN| derivative at $x_5=5$ gets replaced with the\r\n% average slope $d_5 =\\left(\\delta_4 +\\delta_5 \\right)\/2=0.5$.\r\n% \r\n% But we are still not ready!\r\n% \r\n% When should we switch from Akima's formula to the averaging formula? What \r\n% if the slopes are _almost_ equal?\r\n% \r\n% Let's change the above example such that $\\delta_5 =\\varepsilon$ and\r\n% $\\delta_6 =-\\varepsilon$ by setting $v_6 =1+\\varepsilon$ , where\r\n% $\\varepsilon =2^{-52}$ represents\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/eps.html machine epsilon>. We\r\n% can now compare Akima's full formula with the averaging formula at $d_5$\r\n% for two data sets that differ only by $\\varepsilon$:\r\n\r\n% v2eps: new data set which differs only by an eps in v2(6)\r\n    v2eps = v2; v2eps(6) = v2eps(6) + eps;\r\n\r\n    plot(x2,v2,'ko','LineWidth',2,'MarkerSize',10,'DisplayName','Input data v2')\r\n    hold on\r\n    plot(xq2,akima(x2,v2,xq2),'Color',[1 2\/3 0],'LineWidth',2,...\r\n        'DisplayName','v2(6) = 1 uses averaging formula')\r\n    plot(xq2,akima(x2,v2eps,xq2),'-.','Color',[1 2\/3 0],'LineWidth',2,...\r\n        'DisplayName','v2(6) = 1+eps uses Akima''s formula')\r\n    hold off\r\n    legend('Location','SE')\r\n    title('Akima interpolant for (almost) equal side slopes' )\r\n    ylim([-1.2 1.2])\r\n    \r\n%% \r\n% The averaging formula (solid line) returns the derivative $d_5=0.5$ while\r\n% Akima's formula (dashed line) returns $d_5=1$. Therefore, there is an\r\n% unwanted non-negligible difference at around $x_5=5$ between the two\r\n% Akima interpolants, even though the two underlying data sets differ only\r\n% by $\\varepsilon$ at $v_6$.\r\n% \r\n% *Requirement (1)* \r\n%\r\n% We do not want our Akima interpolant to switch to a different \r\n% formula for edge cases.\r\n%\r\n% *Overshoot if data is constant for more than two consecutive nodes*\r\n%\r\n% The previous example was strategically chosen to also reveal another\r\n% property of the Akima interpolant: if the data is constant for more than\r\n% two consecutive nodes, like $v_5=v_6=v_7=1$ in the previous example,\r\n% then the Akima interpolant may produce an overshoot, namely, the\r\n% interpolant's undulation in the $[5~6)$ interval above.\r\n% \r\n% This special type of overshoot is not desirable in many engineering\r\n% applications.\r\n% \r\n% The same concerns hold for the undershoot present in the $[2~3)$ interval \r\n% above.\r\n% \r\n% *Requirement (2)*\r\n%\r\n% We do not want our Akima interpolant to produce overshoot (or undershoot)\r\n% when the data is constant for more than two consecutive nodes.\r\n%\r\n%% Modified Akima formula\r\n% In MATLAB, |'makima'| cubic Hermite interpolation addresses requirements\r\n% (1) and (2) outlined above.\r\n% \r\n% To eliminate overshoot and avoid edge cases of both numerator and\r\n% denominator being equal to 0, we modify Akima's derivative formula by\r\n% tweaking the weights $w_1$ and $w_2$ of the slopes $\\delta_{i-1}$ and\r\n% $\\delta_i$:\r\n% \r\n% $$\r\n% d_i = \\frac{w_1}{w_1+w_2}\\delta_{i-1} + \\frac{w_2}{w_1+w_2}\\delta_i \\,, \r\n% \\quad\\quad w_1 = |\\delta_{i+1} - \\delta_i| + |\\delta_{i+1} + \\delta_i|\/2,\r\n% \\quad\\quad\r\n% w_2 = |\\delta_{i-1} - \\delta_{i-2}| + |\\delta_{i-1} + \\delta_{i-2}| \/ 2.\r\n% $$\r\n%\r\n% _*This modified formula represents the |'makima'| derivative formula used\r\n% in MATLAB:*_\r\n%\r\n% * The addition of the $\\left|\\delta_{i+1} +\\delta_i \\right|\/2$ and\r\n% $\\left|\\delta_{i-1} +\\delta_{i-2} \\right|\/2$ terms forces $d_i =0$ when\r\n% $\\delta_i =\\delta_{i+1} =0$, i.e., $d_i =0$ when $v_i =v_{i+1} =v_{i+2}$.\r\n% Therefore, it eliminates overshoot when the data is constant for more\r\n% than two consecutive nodes.\r\n% * These new terms also address the edge case of equal side slopes\r\n% discussed above, $\\delta_{i-2} =\\delta_{i-1}$ and $\\delta_i\r\n% =\\delta_{i+1}$. However, there is one case which slips through: if the\r\n% data is constant $v_{i-2} = v_{i-1} = v_i = v_{i+1} = v_{i+2}$, then the\r\n% four slopes are all equal to zero and we get $d_i = \\mathrm{NaN}$. For\r\n% this special case of constant data, we set $d_i =0$.\r\n%\r\n% Let's try the |'makima'| formula on the above overshoot example:\r\n\r\n    compareCubicPlots(x2,v2,xq2,true,'SE')\r\n    ylim([-1.2 1.2])\r\n    \r\n%% \r\n% Indeed, |'makima'| does not produce an overshoot if the data is constant\r\n% for more than two nodes ($v_5=v_6=v_7=1$ above).\r\n% \r\n% But what does this mean for the undulations we saw in our first example?\r\n\r\n    compareCubicPlots(x1,v1,xq1,true,'NE')\r\n    \r\n%% \r\n% Notice that |'makima'| closely follows the result obtained with Akima's\r\n% formula. In fact, the results are so similar that it is hard to tell them\r\n% apart on the plot.\r\n% \r\n% Therefore, |'makima'| still preserves Akima's desirable properties of\r\n% being a nice middle ground between |'spline'| and |'pchip'| in terms of\r\n% the resulting undulations.\r\n%\r\n%% Generalization to n-D grids\r\n%\r\n% But we are still not done!\r\n% \r\n% Akima's formula and our modified |'makima'| formula have another\r\n% desirable property: they generalize to higher dimensional n-D gridded\r\n% data. Akima noticed this property in his 1974\r\n% <https:\/\/dl.acm.org\/citation.cfm?id=360779 paper>.\r\n%\r\n% _H. Akima, \"A Method of Bivariate Interpolation and Smooth Surface\r\n% Fitting Based on Local Procedures\", Communications of the ACM, v. 17-1,\r\n% p.18-20, 1974._\r\n%\r\n% For example, to interpolate on a 2-D grid $\\left(x,y\\right)$, we first\r\n% apply the |'makima'| derivative formula separately along $x$ and $y$ to\r\n% obtain two directional derivatives for each grid node. Then, we also\r\n% compute the cross-derivative along $xy$ for each grid node. The\r\n% cross-derivative formula first computes the 2-D divided differences\r\n% corresponding to the 2-D grid data and applies the |'makima'| derivative\r\n% formula along $x$ on these differences; then, it takes the result and\r\n% applies the derivative formula along $y$. The derivatives and\r\n% cross-derivatives are then plugged in as coefficients of a two-variable\r\n% cubic Hermite polynomial representing the 2-D interpolant.\r\n% \r\n% The same procedure applies to higher dimensional n-D grids with $n\\ge 2$\r\n% and requires computing cross-derivatives along all possible\r\n% cross-directions. Therefore, |'makima'| is supported not only in\r\n% |interp1|, but also in |interp2|, |interp3|, |interpn|, and\r\n% |griddedInterpolant|.\r\n% \r\n% Here's how 2-D |'makima'| interpolation compares with 2-D |'spline'|\r\n% interpolation on gridded data generated with the |peaks| function:\r\n\r\n    [X1,Y1,V1] = peaks(5);\r\n    [Xq1,Yq1] = meshgrid(-3:0.1:3,-3:0.1:3);\r\n\r\n    Vqs1 = interp2(X1,Y1,V1,Xq1,Yq1,'spline');\r\n    surf(Xq1,Yq1,Vqs1)\r\n    axis tight\r\n    title('2-D ''spline''')\r\n    snapnow\r\n\r\n    Vqm1 = interp2(X1,Y1,V1,Xq1,Yq1,'makima');\r\n    surf(Xq1,Yq1,Vqm1)\r\n    axis tight\r\n    title('2-D ''makima''')\r\n    \r\n%% \r\n% Notice the smaller undulations (or wiggles) generated by |'makima'|.\r\n% \r\n% Finally, let's try |'makima'| on an example with a few 2-D peaks where\r\n% the data has sharp edges and steps:\r\n\r\n    V2 = zeros(10);\r\n    V2(2:5,2:5) = 3\/7; % First step\r\n    V2(6:7,6:7) = [1\/4 1\/5; 1\/5 1\/4]; % Middle step\r\n    V2(8:9,8:9) = 1\/2; % Last step\r\n    V2 = flip(V2,2);\r\n    [Xq2,Yq2] = meshgrid(1:0.2:10,1:0.2:10);\r\n    \r\n    Vqs2 = interp2(V2,Xq2,Yq2,'spline');\r\n    surf(Xq2,Yq2,Vqs2)\r\n    axis tight\r\n    title('2-D ''spline''')\r\n    snapnow\r\n    \r\n    Vqm2 = interp2(V2,Xq2,Yq2,'makima');\r\n    surf(Xq2,Yq2,Vqm2)\r\n    axis tight\r\n    title('2-D ''makima''')\r\n    \r\n%% \r\n% The |'makima'| result has smaller undulations (or wiggles) than\r\n% |'spline'.|\r\n%\r\n%% Summary\r\n% We conclude with a summary of properties of interest for the three cubic\r\n% Hermite interpolation methods supported in MATLAB:\r\n\r\n    summary = table( ...\r\n        [\"C2\"; \"Uses all points\"; \"Not-a-knot condition\"; \"Yes\"], ...\r\n        [\"C1\"; \"Uses 4 points\"; \"One-sided slope\"; \"No\"], ...\r\n        [\"C1\"; \"Uses 5 points\"; \"Quadratic extrapolation\"; \"Yes\"], ...\r\n        'VariableNames', ...\r\n        [\"spline\" \"pchip\" \"makima\"], ...\r\n        'RowNames', ...\r\n        [\"Continuity\"; \"Locality\"; \"End points\"; \"Supports n-D\"]);\r\n    disp(summary)\r\n\r\n \r\n%% |code|\r\n\r\n%% |compareCubicPlots|\r\nfunction compareCubicPlots(x,v,xq,showMakima,legendLocation)\r\n% Plot 'pchip', 'spline', Akima, and 'makima' interpolation results\r\n    vqp = pchip(x,v,xq);      % same as vq = interp1(x,v,xq,'pchip')\r\n    vqs = spline(x,v,xq);     % same as vq = interp1(x,v,xq,'spline')\r\n    vqa = akima(x,v,xq);\r\n    if showMakima\r\n        vqm = makima(x,v,xq); % same as vq = interp1(x,v,xq,'makima')\r\n    end\r\n\r\n    plot(x,v,'ko','LineWidth',2,'MarkerSize',10,'DisplayName','Input data')\r\n    hold on\r\n    plot(xq,vqp,'LineWidth',4,'DisplayName','''pchip''')\r\n    plot(xq,vqs,'LineWidth',2,'DisplayName','''spline''')\r\n    plot(xq,vqa,'LineWidth',2,'DisplayName','Akima''s formula')\r\n    if showMakima\r\n        plot(xq,vqm,'REPLACE_WITH_DASH_DASH','Color',[0 3\/4 0],'LineWidth',2, ...\r\n            'DisplayName','''makima''')\r\n    end\r\n    hold off\r\n    xticks(x(1)-1:x(end)+1)\r\n    legend('Location',legendLocation)\r\n    title('Cubic Hermite interpolation')\r\nend\r\n\r\n%% |akima|\r\nfunction vq = akima(x,v,xq)\r\n%AKIMA   Akima piecewise cubic Hermite interpolation.\r\n%\r\n%   vq = AKIMA(x,v,xq) interpolates the data (x,v) at the query points xq\r\n%   using Akima's cubic Hermite interpolation formula.\r\n%\r\n%   References:\r\n%\r\n%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting\r\n%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.\r\n%\r\n%   AKIMA vs. PCHIP vs. SPLINE:\r\n%\r\n%     - Akima's cubic formula is a middle ground between SPLINE and PCHIP:\r\n%       It has lower-amplitude wiggles than SPLINE, but is not as agressive\r\n%       at reducing the wiggles as PCHIP.\r\n%     - Akima's cubic formula and SPLINE generalize to n-D grids.\r\n%\r\n%   Example: Compare AKIMA with MAKIMA, SPLINE, and PCHIP\r\n%\r\n%       x = [1 2 3 4 5 5.5 7 8 9 9.5 10];\r\n%       v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];\r\n%       xq = 0.75:0.05:10.25;\r\n%       vqs = spline(x,v,xq);\r\n%       vqp = pchip(x,v,xq);\r\n%       vqa = akima(x,v,xq);\r\n%       vqm = makima(x,v,xq);\r\n%\r\n%       figure\r\n%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on\r\n%       plot(xq,vqp,'LineWidth',4)\r\n%       plot(xq,vqs,xq,vqa,xq,vqm,'REPLACE_WITH_DASH_DASH','LineWidth',2)\r\n%       legend('Data','''pchip''','''spline''','Akima''s formula','''makima''')\r\n%       title('Cubic Hermite interpolation in MATLAB')\r\n%\r\n%   See also MAKIMA, SPLINE, PCHIP.\r\n\r\n%   Copyright 2019 The MathWorks, Inc.\r\n%   mailto: cosmin.ionita@mathworks.com\r\n\r\n    assert(isvector(x) && isvector(v) && (numel(x) == numel(v)))\r\n    assert(numel(x) > 2)\r\n    x = x(:).'; v = v(:).'; % Same shapes as in pchip.m and spline.m\r\n\r\n    % Compute Akima slopes\r\n    h = diff(x);\r\n    delta = diff(v).\/h;\r\n    slopes = akimaSlopes(delta);\r\n\r\n    % Evaluate the piecewise cubic Hermite interpolant\r\n    vq = ppval(pwch(x,v,slopes,h,delta),xq);\r\nend\r\n\r\n%% |akimaSlopes|\r\nfunction s = akimaSlopes(delta)\r\n% Derivative values for Akima cubic Hermite interpolation\r\n\r\n% Akima's derivative estimate at grid node x(i) requires the four finite\r\n% differences corresponding to the five grid nodes x(i-2:i+2).\r\n%\r\n% For boundary grid nodes x(1:2) and x(n-1:n), append finite differences\r\n% which would correspond to x(-1:0) and x(n+1:n+2) by using the following\r\n% uncentered difference formula correspondin to quadratic extrapolation\r\n% using the quadratic polynomial defined by data at x(1:3)\r\n% (section 2.3 in Akima's paper):\r\n    n = numel(delta) + 1; % number of grid nodes x\r\n    delta_0  = 2*delta(1)   - delta(2);\r\n    delta_m1 = 2*delta_0    - delta(1);\r\n    delta_n  = 2*delta(n-1) - delta(n-2);\r\n    delta_n1 = 2*delta_n    - delta(n-1);\r\n    delta = [delta_m1 delta_0 delta delta_n delta_n1];\r\n\r\n% Akima's derivative estimate formula (equation (1) in the paper):\r\n%\r\n%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting\r\n%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.\r\n%\r\n% s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i))\r\n%      \/ (|d(i+1)-d(i)|          + |d(i-1)-d(i-2)|)\r\n    weights = abs(diff(delta));\r\n    weights1 = weights(1:n);   % |d(i-1)-d(i-2)|\r\n    weights2 = weights(3:end); % |d(i+1)-d(i)| \r\n    delta1 = delta(2:n+1);     % d(i-1)\r\n    delta2 = delta(3:n+2);     % d(i)\r\n\r\n    weights12 = weights1 + weights2;\r\n    s = (weights2.\/weights12) .* delta1 + (weights1.\/weights12) .* delta2;\r\n\r\n% To avoid 0\/0, Akima proposed to average the divided differences d(i-1)\r\n% and d(i) for the edge case of d(i-2) = d(i-1) and d(i) = d(i+1):\r\n    ind = weights1 == 0 & weights2 == 0;\r\n    s(ind) = (delta1(ind) + delta2(ind)) \/ 2;\r\nend\r\n\r\n%% |makima|\r\nfunction vq = makima(x,v,xq)\r\n%MAKIMA   Modified Akima piecewise cubic Hermite interpolation.\r\n%\r\n%       We modify Akima's formula to eliminate overshoot and undershoot\r\n%       when the data is constant for more than two consecutive nodes.\r\n%\r\n%   vq = MAKIMA(x,v,xq) interpolates the data (x,v) at the query points xq\r\n%   using a modified Akima cubic Hermite interpolation formula.\r\n%\r\n%   References:\r\n%\r\n%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting\r\n%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.\r\n%\r\n%   MAKIMA vs. PCHIP vs. SPLINE:\r\n%\r\n%     - MAKIMA is a middle ground between SPLINE and PCHIP:\r\n%       It has lower-amplitude wiggles than SPLINE, but is not as agressive\r\n%       at reducing the wiggles as PCHIP.\r\n%     - MAKIMA and SPLINE generalize to n-D grids.\r\n%\r\n%   Example: No overshoot and undershoot with MAKIMA when the data is\r\n%            constant for more than two consecutive nodes\r\n%\r\n%       x = 1:7;\r\n%       v = [-1 -1 -1 0 1 1 1];\r\n%       xq = 0.75:0.05:7.25;\r\n%       vqs = spline(x,v,xq);\r\n%       vqp = pchip(x,v,xq);\r\n%       vqa = akima(x,v,xq);\r\n%       vqm = makima(x,v,xq);\r\n%\r\n%       figure\r\n%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on\r\n%       plot(xq,vqp,'LineWidth',4)\r\n%       plot(xq,vqs,xq,vqa,xq,vqm,'LineWidth',2)\r\n%       legend('Data','''pchip''','''spline''','Akima''s formula',...\r\n%           '''makima''','Location','SE')\r\n%       title('''makima'' has no overshoot in x(1:3) and x(5:7)')\r\n%\r\n%   Example: Compare MAKIMA with AKIMA, SPLINE, and PCHIP\r\n%\r\n%       x = [1 2 3 4 5 5.5 7 8 9 9.5 10];\r\n%       v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];\r\n%       xq = 0.75:0.05:10.25;\r\n%       vqs = spline(x,v,xq);\r\n%       vqp = pchip(x,v,xq);\r\n%       vqa = akima(x,v,xq);\r\n%       vqm = makima(x,v,xq);\r\n%\r\n%       figure\r\n%       plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on\r\n%       plot(xq,vqp,'LineWidth',4)\r\n%       plot(xq,vqs,xq,vqa,xq,vqm,'REPLACE_WITH_DASH_DASH','LineWidth',2)\r\n%       legend('Data','''pchip''','''spline''','Akima''s formula','''makima''')\r\n%       title('Cubic Hermite interpolation in MATLAB')\r\n%\r\n%   See also AKIMA, SPLINE, PCHIP.\r\n\r\n%   Copyright 2019 The MathWorks, Inc.\r\n%   mailto: cosmin.ionita@mathworks.com\r\n\r\n    assert(isvector(x) && isvector(v) && (numel(x) == numel(v)))\r\n    assert(numel(x) > 2)\r\n    x = x(:).'; v = v(:).'; % Same shapes as in pchip.m and spline.m\r\n\r\n% Compute modified Akima slopes\r\n    h = diff(x);\r\n    delta = diff(v).\/h;\r\n    slopes = makimaSlopes(delta);\r\n\r\n% Evaluate the piecewise cubic Hermite interpolant\r\n    vq = ppval(pwch(x,v,slopes,h,delta),xq);\r\n   \r\nend\r\n\r\n%% |makimaSlopes|\r\nfunction s = makimaSlopes(delta)\r\n% Derivative values for modified Akima cubic Hermite interpolation\r\n\r\n% Akima's derivative estimate at grid node x(i) requires the four finite\r\n% differences corresponding to the five grid nodes x(i-2:i+2).\r\n%\r\n% For boundary grid nodes x(1:2) and x(n-1:n), append finite differences\r\n% which would correspond to x(-1:0) and x(n+1:n+2) by using the following\r\n% uncentered difference formula correspondin to quadratic extrapolation\r\n% using the quadratic polynomial defined by data at x(1:3)\r\n% (section 2.3 in Akima's paper):\r\n    n = numel(delta) + 1; % number of grid nodes x\r\n    delta_0  = 2*delta(1)   - delta(2);\r\n    delta_m1 = 2*delta_0    - delta(1);\r\n    delta_n  = 2*delta(n-1) - delta(n-2);\r\n    delta_n1 = 2*delta_n    - delta(n-1);\r\n    delta = [delta_m1 delta_0 delta delta_n delta_n1];\r\n\r\n% Akima's derivative estimate formula (equation (1) in the paper):\r\n%\r\n%       H. Akima, \"A New Method of Interpolation and Smooth Curve Fitting\r\n%       Based on Local Procedures\", JACM, v. 17-4, p.589-602, 1970.\r\n%\r\n% s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i))\r\n%      \/ (|d(i+1)-d(i)|          + |d(i-1)-d(i-2)|)\r\n%\r\n% To eliminate overshoot and undershoot when the data is constant for more\r\n% than two consecutive nodes, in MATLAB's 'makima' we modify Akima's\r\n% formula by adding an additional averaging term in the weights.\r\n% s(i) = ( (|d(i+1)-d(i)|   + |d(i+1)+d(i)|\/2  ) * d(i-1) +\r\n%          (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|\/2) * d(i)  )\r\n%      \/ ( (|d(i+1)-d(i)|   + |d(i+1)+d(i)|\/2  ) +\r\n%          (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|\/2)\r\n    weights = abs(diff(delta)) + abs((delta(1:end-1)+delta(2:end))\/2);\r\n\r\n    weights1 = weights(1:n);   % |d(i-1)-d(i-2)|\r\n    weights2 = weights(3:end); % |d(i+1)-d(i)| \r\n    delta1 = delta(2:n+1);     % d(i-1)\r\n    delta2 = delta(3:n+2);     % d(i)\r\n\r\n    weights12 = weights1 + weights2;\r\n    s = (weights2.\/weights12) .* delta1 + (weights1.\/weights12) .* delta2;\r\n\r\n% If the data is constant for more than four consecutive nodes, then the\r\n% denominator is zero and the formula produces an unwanted NaN result.\r\n% Replace this NaN with 0.\r\n    s(weights12 == 0) = 0;\r\nend\r\n##### SOURCE END ##### e5dc3d9c091f41a6acd2477ff3ef9b53\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/makima_image.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction--><p><i>This post is by my colleague <b>Cosmin Ionita<\/b>.<\/i>... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2019\/04\/29\/makima-piecewise-cubic-interpolation\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":4859,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,12,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/4707"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=4707"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/4707\/revisions"}],"predecessor-version":[{"id":4861,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/4707\/revisions\/4861"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/4859"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=4707"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=4707"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=4707"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}