📄 csapidem.html
字号:
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Cubic spline interpolation</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="csapidem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit csapidem">Open csapidem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo csapidem">Run in the Command Window</a></div> </div> <div class="content"> <h1>Cubic spline interpolation</h1> <introduction> <p>Illustration of the use of CSAPI and CSAPE.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">CSAPI</a></li> <li><a href="#2">Two data points</a></li> <li><a href="#3">Three data points</a></li> <li><a href="#4">Four or more points</a></li> <li><a href="#5">How to check CSAPI</a></li> <li><a href="#6">The truncated power</a></li> <li><a href="#7">A particular truncated power</a></li> <li><a href="#8">The particular truncated power interpolated</a></li> <li><a href="#9">The error</a></li> <li><a href="#10">A truncated power that cannot be reproduced</a></li> <li><a href="#12">PPFORM instead of values</a></li> <li><a href="#13">Differentiating the interpolant</a></li> <li><a href="#14">Error</a></li> <li><a href="#15">Second derivative</a></li> <li><a href="#16">Integral</a></li> <li><a href="#17">CSAPE</a></li> <li><a href="#18">csapi(x,y) vs csape(x,y)</a></li> <li><a href="#19">Other end conditions: the `natural' spline interpolant</a></li> <li><a href="#20">The `natural' spline interpolant: an example</a></li> <li><a href="#21">Other end conditions: prescribing second derivatives</a></li> <li><a href="#22">Other end conditions: prescribing slopes</a></li> <li><a href="#23">Other end conditions: mixed end conditions</a></li> <li><a href="#24">Other end conditions: periodic conditions</a></li> <li><a href="#25">End conditions not explicitly covered by CSAPI or CSAPE</a></li> </ul> </div> <h2>CSAPI<a name="1"></a></h2> <p>The spline toolbox command</p><pre> values = csapi(x,y,xx)</pre><p>returns the values at XX of the cubic spline interpolant to the given data (X,Y), using the not-a-knot end condition. This interpolant is a piecewise cubic with break sequence X whose cubic pieces join together to form a function with two continuous derivatives. The "not-a-knot" end condition means that, at the first and last interior break, even the third derivative is continuous (up to roundoff). </p> <h2>Two data points<a name="2"></a></h2> <p>Just two data points result in a straight line interpolant:</p><pre class="codeinput">xx = linspace(0,6,121); x = [0 1]; y = [2 0];plot(xx,csapi(x,y,xx),<span class="string">'k-'</span>,x,y,<span class="string">'ro'</span>), grid <span class="string">off</span>title(<span class="string">'Interpolant to two points'</span>)</pre><img vspace="5" hspace="5" src="csapidem_01.png"> <h2>Three data points<a name="3"></a></h2> <p>Three data points give a parabola:</p><pre class="codeinput">x = [2 3 5]; y = [1 0 4];plot(xx,csapi(x,y,xx),<span class="string">'k-'</span>,x,y,<span class="string">'ro'</span>), grid <span class="string">off</span>title(<span class="string">'Interpolant to three points'</span>)</pre><img vspace="5" hspace="5" src="csapidem_02.png"> <h2>Four or more points<a name="4"></a></h2> <p>Four or more data points give in general a cubic spline.</p><pre class="codeinput">x = [1 1.5 2 4.1 5]; y = [1 -1 1 -1 1];plot(xx,csapi(x,y,xx),<span class="string">'k-'</span>,x,y,<span class="string">'ro'</span>), grid <span class="string">off</span>title(<span class="string">'Cubic spline interpolant to five points'</span>)</pre><img vspace="5" hspace="5" src="csapidem_03.png"> <h2>How to check CSAPI<a name="5"></a></h2> <p>These look like nice interpolants, but how do we check that CSAPI performs as advertised?</p> <p>We already saw that we interpolate, for we always plotted the data points and the interpolant went right through those points.</p> <p>But, to be sure that we get a cubic spline, it is best to start with data from a cubic spline of the expected sort and check whether CSAPI gives us back that cubic spline. </p> <h2>The truncated power<a name="6"></a></h2> <p>One such is the truncated third power, i.e., the function x |-> subplus(x - xi)^3 , with XI one of the breaks and SUBPLUS the t r u n c a t i o n f u n c t i o n, provided by the spline toolbox command SUBPLUS: </p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span>help <span class="string">subplus</span></pre><pre class="codeoutput"> SUBPLUS Positive part. x , if x>=0 y = subplus(x) := (x)_{+} = , 0 , if x<=0 returns the positive part of X. Used for computing truncated powers.</pre><img vspace="5" hspace="5" src="csapidem_04.png"> <h2>A particular truncated power<a name="7"></a></h2> <p>For the particular choice XI = 2 , we get the function plotted above. As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right of 2. </p><pre class="codeinput">plot(xx, subplus(xx-2).^3,<span class="string">'y'</span>,<span class="string">'linewidth'</span>,3), grid <span class="string">off</span>axis([0,6,-10,70])</pre><img vspace="5" hspace="5" src="csapidem_05.png"> <h2>The particular truncated power interpolated<a name="8"></a></h2> <p>Now we interpolate this particular cubic spline at the data sites 0:6, and plot the interpolant on top of the spline, in black.</p><pre class="codeinput">x = 0:6; y = subplus(x-2).^3;values = csapi(x,y,xx);hold <span class="string">on</span>, plot(xx,values,<span class="string">'k'</span>,x,y,<span class="string">'ro'</span>), hold <span class="string">off</span>title(<span class="string">'Interpolant to subplus(x-2)^3'</span>)</pre><img vspace="5" hspace="5" src="csapidem_06.png"> <h2>The error<a name="9"></a></h2> <p>When comparing two functions, it is usually much more informative to plot their difference. To put the size of their difference into context, we also compute the maximum data value. This shows the error to be no worse than the inevitable round-off. </p><pre class="codeinput">plot(xx,values-subplus(xx-2).^3)title(<span class="string">'Error in cubic spline interpolation to subplus(x-2)^3'</span>)max_y=max(abs(y))</pre><pre class="codeoutput">max_y = 64</pre><img vspace="5" hspace="5" src="csapidem_07.png"> <h2>A truncated power that cannot be reproduced<a name="10"></a></h2> <p>As a further test, we interpolate a truncated power whose CSAPI interpolant at the sites 0:6 cannot coincide with it. For example, the first interior break of the interpolating spline is not really a knot since the interpolant has three continuous derivatives there. Hence we should not be able to reproduce the truncated power centered at that site. We try </p><pre class="codeinput">values = csapi(x,subplus(x-1).^3,xx);plot(xx,values-subplus(xx-1).^3)title(<span class="string">'Error in not-a-knot interpolant to subplus(x-1)^3'</span>)xlabel([<span class="string">'since 1 is a first interior knot, it is not active for this'</span>,<span class="keyword">...</span> <span class="string">' interpolant'</span>])</pre><img vspace="5" hspace="5" src="csapidem_08.png"> <p>The difference is as large as .18, but decays rapidly as we move away from 1. This illustrates that cubic spline interpolation is essentially local. </p> <h2>PPFORM instead of values<a name="12"></a></h2> <p>It is possible to retain the interpolating cubic spline in a form suitable for subsequent evaluation, or for calculating its derivatives, or for other manipulations. This is done by calling CSAPI in the form </p><pre> pp = csapi(x,y)</pre><p>which returns, in PP, the ppform of the interpolant. This form can be evaluated at some points XX by</p><pre> values = fnval(pp,xx)</pre><p>It can be differentiated by</p><pre> dpp = fnder(pp)</pre><p>or integrated by</p><pre> ipp = fnint(pp)</pre><p>which return, in DPP or IPP, the ppform of the derivative or the primitive, respectively.</p> <h2>Differentiating the interpolant<a name="13"></a></h2> <p>For example, we plot the derivative of this truncated power, i.e., x |-> 3*subplus(x-2)^2 (again in yellow) and, on top of it, the derivative of our interpolant (again in black). </p><pre class="codeinput">pp = csapi(x,subplus(x-2).^3); dpp = fnder(pp);plot(xx,fnval(dpp,xx),<span class="string">'y'</span>,<span class="string">'linewidth'</span>,3), hold <span class="string">on</span>plot(xx,3*subplus(xx-2).^2,<span class="string">'k'</span>), hold <span class="string">off</span>, grid <span class="string">off</span>title(<span class="string">'Derivative of interpolant to subplus(x-2)^3'</span>)</pre><img vspace="5" hspace="5" src="csapidem_09.png"> <h2>Error<a name="14"></a></h2> <p>Again, the more informative comparison is to plot their difference, and this is again no bigger than round-off.</p><pre class="codeinput">plot(xx,fnval(dpp,xx)-3*subplus(xx-2).^2)title(<span class="string">'(derivative of interpolant to subplus(x-2)^3 ) - 3*subplus(x-2)^2'</span>)</pre><img vspace="5" hspace="5" src="csapidem_10.png"> <h2>Second derivative<a name="15"></a></h2> <p>The second derivative of the interpolant should be 6*subplus(.-2). We try it, plotting above the difference between this function and the second derivative of the interpolant to subplus(.-2)^3. Now there are jumps, but they are still within roundoff. </p><pre class="codeinput">ddpp = fnder(dpp); plot(xx,fnval(ddpp,xx)-6*subplus(xx-2))title(<span class="string">'(second derivative of interpolant to subplus(x-2)^3 ) - 6*subplus(x-2)'</span>)</pre><img vspace="5" hspace="5" src="csapidem_11.png"> <h2>Integral<a name="16"></a></h2> <p>The integral of subplus(.-2)^3 is subplus(.-2)^4/4. We plot the difference between it and the integral of the interpolant to subplus(.-2)^3: </p><pre class="codeinput">ipp = fnint(pp);plot(xx,fnval(ipp,xx)-subplus(xx-2).^4/4)title(<span class="string">'(integral of interpolant to subplus(x-2)^3 ) - subplus(x-2)^4/4'</span>)</pre><img vspace="5" hspace="5" src="csapidem_12.png"> <h2>CSAPE<a name="17"></a></h2> <p>The command CSAPE also provides a cubic spline interpolant to given data, but permits various other end conditions. It does not directly return values of that interpolant, but only its ppform. Its simplest version, </p><pre> pp = csape(x,y)</pre><p>uses the Lagrange end conditions, which are better at times than the not-a-knot condition used by CSAPI.</p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapidem_13.png"> <h2>csapi(x,y) vs csape(x,y)<a name="18"></a></h2> <p>For example, consider again interpolation to the truncated power which CSAPI fails to reproduce. We plot here (in black), the error of the not-a-knot interpolant, along with the error (in red) of the interpolant obtained from CSAPE; -- not much difference between the two in this case. </p><pre class="codeinput">exact = subplus(xx-1).^3;plot(xx, fnval( csapi(x,subplus(x-1).^3), xx ) - exact,<span class="string">'k'</span>), hold <span class="string">on</span>plot(xx, fnval( csape(x,subplus(x-1).^3), xx ) - exact,<span class="string">'r'</span>)title(<span class="string">'Error in not-a-knot (black) vs. Lagrange (red)'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapidem_14.png"> <h2>Other end conditions: the `natural' spline interpolant<a name="19"></a></h2> <p>The command CSAPE only provides the interpolating cubic spline in ppform, but for several different end conditions. For example, the call </p><pre> pp = csape(x,y,'variational')</pre><p>uses the so-called `natural' end conditions. This means that the second derivative is zero at the two extreme breaks.</p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapidem_15.png"> <h2>The `natural' spline interpolant: an example<a name="20"></a></h2> <p>We apply `natural' cubic spline interpolation to the truncated power we were able to reproduce, and plot the error. We now get a large error near the right end, due to the fact that we insisted on a zero second derivative there. For variety, we get the `natural' spline interpolant here by the command csape(x,y,'second') which uses the default value 0 for the second derivative at the extreme data sites. </p><pre class="codeinput">pp = csape(x,subplus(x-2).^3,<span class="string">'second'</span>);plot(xx, fnval(pp,xx) - subplus(xx-2).^3 )title(<span class="string">'Error in ''natural'' spline interpolation to subplus(x-2)^3'</span>)xlabel(<span class="string">' note the large error near the right end'</span>)</pre><img vspace="5" hspace="5" src="csapidem_16.png"> <h2>Other end conditions: prescribing second derivatives<a name="21"></a></h2> <p>When we use explicitly the correct second derivatives, we get a small error, as the above plot shows. The command CSAPE(X,[VALS(1),Y VALS(2)],'second') specifies that second derivatives are to be matched, with VALS(1) the second derivative at the left endpoint, and VALS(2) the second derivative at the right. We compute these values directly from the second derivative of the truncated power subplus(.-2)^3. </p><pre class="codeinput">pp = csape(x,[6*subplus(x(1)-2),subplus(x-2).^3,6*subplus(x(end)-2)], <span class="string">'second'</span>);plot(xx, fnval( pp, xx ) - subplus(xx-2).^3,<span class="string">'r'</span>)title(<span class="string">'Error in spline interpolation to subplus(x-1)^3 when matching f'''' at ends '</span>)</pre><img vspace="5" hspace="5" src="csapidem_17.png"> <h2>Other end conditions: prescribing slopes<a name="22"></a></h2> <p>CSAPE also permits specification of endpoint s l o p e s . This is the c l a m p e d spline (or, c o m p l e t e cubic spline interpolant). The statement </p><pre> pp = csape(x,[sl,y,sr],'clamped')</pre><p>supplies the cubic spline interpolant to the data X, Y that also has slope SL at the leftmost data site and slope SR at the rightmost data site. </p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapidem_18.png"> <h2>Other end conditions: mixed end conditions<a name="23"></a></h2> <p>It is even possible to mix these conditions. For example, our much exercised truncated power function f(x) = subplus(x-1)^3 has slope 0 at x = 0 and second derivative 30 at x = 6 (our last data site). </p> <p>We therefore expect no error in the following interpolant:</p><pre class="codeinput">pp = csape(x, [0,subplus(x-1).^3,30], [1 2] );plot(xx, fnval(pp,xx) - subplus(xx-1).^3 )title(<span class="string">'Error in spline interpolation to subplus(x-1)^3 when matching ...'</span>)xlabel(<span class="string">' ... slope at left end and curvature at right.'</span>)</pre><img vspace="5" hspace="5" src="csapidem_19.png"> <h2>Other end conditions: periodic conditions<a name="24"></a></h2> <p>It is also possible to prescribe p e r i o d i c end conditions. For example, the sine function is 2*pi periodic and has the values [0 -1 0 1 0] at the sites (pi/2)*(-2:2) . The difference, plotted above, between the sine function and its periodic cubic spline interpolant at these sites, is only 2 percent. Not bad. </p><pre class="codeinput">x = (pi/2)*(-2:2);y = [0 -1 0 1 0];pp = csape(x,y, <span class="string">'periodic'</span> );xx = linspace(-pi,pi,201);plot(xx, sin(xx) - fnval(pp,xx), <span class="string">'x'</span>)title(<span class="string">'Error in periodic cubic spline interpolation to sin(x)'</span>)</pre><img vspace="5" hspace="5" src="csapidem_20.png"> <h2>End conditions not explicitly covered by CSAPI or CSAPE<a name="25"></a></h2> <p>Any end condition not covered explicitly by CSAPI or CSAPE can be handled by constructing the interpolant with the CSAPE default side conditions, and then adding to it an appropriate scalar multiple of an interpolant to zero values and some side conditions. If there are two `nonstandard' side conditions to be satisfied, one may have to solve a 2x2 linear system first. </p> <p>For example, suppose that you want to enforce the condition</p><pre> lambda(s) := a (Ds)(e) + b (D^2 s)(e) = c</pre><p>on the cubic spline interpolant s to the following data (taken from a a quartic polynomial that happens to satisfy this specific side condition): </p><pre class="codeinput">cla, title(<span class="string">''</span>), axis <span class="string">off</span>q = inline(<span class="string">'x.*(-1 + x.*(-1+x.*x/5))'</span>);x = 0:.25:3; y = q(x);e = x(1); a = 2; b = -3; c = 4;</pre><img vspace="5" hspace="5" src="csapidem_21.png"> <p>Then, in addition to the interpolant PP1 with the default end conditions ... pp1 = csape(x,y); ... and first derivative DP1 of its first polynomial piece ... dp1 = fnder( fnbrk(pp1,1) ); ... we also construct the cubic spline interpolant PP0 to zero data, and with a slope of 1 at E, as well as the first derivative DP0 of its first polynomial piece. </p><pre class="codeinput">pp1 = csape(x,y);dp1 = fnder(fnbrk(pp1,1));pp0 = csape(x,[1,zeros(size(y)),0], [1,0] );⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -