📄 spapidem.html

📁 演示matlab曲线拟和与插直的基本方法
💻 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>Splines and real world data</title>      <meta name="generator" content="MATLAB 7.1">      <meta name="date" content="2005-07-27">      <meta name="m-file" content="spapidem">      <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css">   </head>   <body>      <div class="header">         <div class="left"><a href="matlab:edit spapidem">Open spapidem.m in the Editor</a></div>         <div class="right"><a href="matlab:echodemo spapidem">Run in the Command Window</a></div>      </div>      <div class="content">         <h1>Splines and real world data</h1>         <introduction>            <p>This a demonstration of a spline function being fit to real-world data.</p>         </introduction>         <h2>Contents</h2>         <div>            <ul>               <li><a href="#1">Manual knot choice for spline interpolation</a></li>               <li><a href="#11">Automatic knot choice for interpolation</a></li>            </ul>         </div>         <h2>Manual knot choice for spline interpolation<a name="1"></a></h2>         <p>Here are some data which record a property of titanium measured as a function of temperature. We'll use it to illustrate some            issues with spline interpolation.         </p>         <p>The plot of the data shows a rather sharp peak.</p><pre class="codeinput">hold <span class="string">off</span>[xx,yy]=titanium;frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];plot(xx,yy,<span class="string">'x'</span>);axis(frame);</pre><img vspace="5" hspace="5" src="spapidem_01.png"> <p>We pick a few data points from these somewhat rough data, since we want to interpolate. Here is a picture of the data, with            the selected data points marked.         </p><pre class="codeinput">hold <span class="string">on</span>pick=[1 5 11 21 27 29 31 33 35 40 45 49];tau=xx(pick);y=yy(pick);plot(tau,y,<span class="string">'ro'</span>);</pre><img vspace="5" hspace="5" src="spapidem_02.png"> <p>Since a spline of order k with n+k knots has n degrees of freedom, and we have 12 data points, a fit with a fourth order spline            requires 12+4 = 16 knots. Moreover, this knot sequence t must be such that the i-th data site lies in the support of the i-th            B-spline. We achieve this by using the data sites as knots, but adding two simple knots at either end.         </p><pre> dl = tau(2)-tau(1); dr = tau(end)-tau(end-1); t = [tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];</pre><p>We use this knot sequence to construct an interpolating cubic spline:</p><pre> sp = spapi(t, tau, y);</pre><pre class="codeinput">dl=tau(2)-tau(1);dr=tau(end)-tau(end-1);t=[tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];   <span class="comment">% construct the knot sequence</span>hold <span class="string">on</span>axis(frame+[-2*dl 2*dr 0 0])plot(t,repmat(frame(3)+.03,size(t)),<span class="string">'kx'</span>)hold <span class="string">off</span>sp=spapi(t,tau,y);                         <span class="comment">% This constructs the spline.</span></pre><img vspace="5" hspace="5" src="spapidem_03.png"> <p>Now, for the plot. Since we do not care about the part of the spline outside the data interval, we restrict the plot to that            interval:         </p><pre>hold on,  fnplt(sp,[tau(1) tau(end)], 'k'),  hold off</pre><pre class="codeinput">plot(xx,yy,<span class="string">'x'</span>,tau,y,<span class="string">'ro'</span>), axis(frame), hold <span class="string">on</span><span class="comment">% Now, for the plot:</span>fnplt(sp,[tau(1) tau(end)], <span class="string">'k'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_04.png"> <p>A closer look at the left part of the spline fit shows some undulations.</p><pre class="codeinput">xxx = linspace(tau(1),tau(5),41);plot(xxx, fnval(sp, xxx), <span class="string">'k'</span>, tau, y, <span class="string">'ro'</span>);axis([tau(1) tau(5) 0.6 1.2]);</pre><img vspace="5" hspace="5" src="spapidem_05.png"> <p>The unreasonable bump in the first interval stems from the fact that our spline goes smoothly to zero at its first knot.</p>         <p>Here is a picture of the entire spline, along with its knot sequence.</p><pre class="codeinput">axis([tau(1) tau(5) 0.6 1.2]);fnplt(sp,<span class="string">'k'</span>);hold <span class="string">on</span>, plot(t,repmat(.1,size(t)),<span class="string">'kx'</span>), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_06.png"> <p>Here are again the data points as well.</p><pre class="codeinput">hold <span class="string">on</span>plot(tau, y, <span class="string">'ro'</span>);hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_07.png"> <p>Here is a simple way to enforce a more reasonable boundary behavior. We add two more data points outside the given data interval            and choose as our data there the values of the straight line through the first two data points.         </p><pre class="codeinput">tt=[tau(1)-[4 3 2 1]*dl tau tau(end)+[1 2 3 4]*dr];xx=[tau(1)-[2 1]*dl tau tau(end)+[1 2]*dr];yy=[y(1)-[2 1]*(y(2)-y(1)) y y(end)+[1 2]*(y(end)-y(end-1))];sp2=spapi(tt,xx,yy); fnplt(sp2,<span class="string">'b'</span>,tau([1 end]))hold <span class="string">on</span>plot(tau,y,<span class="string">'or'</span>, xx([1 2 end-1 end]),yy([1 2 end-1 end]),<span class="string">'ko'</span>);axis(frame+[-2*dl 2*dr 0 0]);hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_08.png"> <p>Here is also the original spline fit, shown in black, to show the reduction of the undulation in the first and last interval.</p><pre class="codeinput">hold <span class="string">on</span>fnplt(sp,<span class="string">'k'</span>,tau([1 end]))hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_09.png"> <p>Finally, here is a closer look at the first four data intervals which shows more clearly the reduction of the undulation near            the left end.         </p><pre class="codeinput">plot(xxx,fnval(sp2,xxx),<span class="string">'b'</span>,tau,y,<span class="string">'ro'</span>,xxx,fnval(sp,xxx),<span class="string">'k'</span>);axis([tau(1) tau(5) .6 2.2]);hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_10.png"> <h2>Automatic knot choice for interpolation<a name="11"></a></h2>         <p>If all this detail turns you off, let the spline toolbox choose the knots for you, by using the spline interpolation command            SPAPI in the form        spapi( order, data_sites, data_values )         </p><pre class="codeinput">autosp = spapi( 4, tau, y);fnplt(autosp,<span class="string">'g'</span>)hold <span class="string">on</span>, plot(tau, y, <span class="string">'or'</span>), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_11.png"> <p>Here is the result of a much better knot choice (plotted above as red x's), obtained by shifting the knot at 842 slightly            to the right and the knot at 985 slightly to the left.         </p><pre class="codeinput">knots = fnbrk(autosp,<span class="string">'knots'</span>);hold <span class="string">on</span>, plot(knots, repmat(.5,size(knots)),<span class="string">'xg'</span>)knots([7 12]) = [851, 971];plot(knots, repmat(.54,size(knots)),<span class="string">'xr'</span>)adjsp = spapi(knots, tau, y);fnplt(adjsp,<span class="string">'r'</span>,2), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_12.png"> <p>Else, simply try the standard cubic spline interpolant, supplied by CSAPI (which amounts to a slightly different choice of            knots):         </p>         <p>The next and last slide shows all five interpolants, for comparison.</p><pre class="codeinput">autocs = csapi(tau, y);fnplt(autocs,<span class="string">'c'</span>)hold <span class="string">on</span>, plot(tau, y, <span class="string">'or'</span>), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_13.png"> <p>With such fast varying data, it is hard to get agreement among all reasonable interpolants, even if each of them is a cubic            spline.         </p><pre class="codeinput">hold <span class="string">on</span>fnplt(sp,<span class="string">'k'</span>,tau([1 end]))  <span class="comment">% black: original</span>fnplt(sp2,<span class="string">'b'</span>,tau([1 end])) <span class="comment">% blue:  with special end conditions</span>fnplt(autosp,<span class="string">'g'</span>)           <span class="comment">% green: automatic knot choice by SPAPI</span>fnplt(autocs,<span class="string">'c'</span>)           <span class="comment">% cyan:  automatic knot choice by CSAPI</span>fnplt(adjsp,<span class="string">'r'</span>,2)          <span class="comment">% red:   knot choice by SPAPI slightly changed</span>hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="spapidem_14.png"> <p class="footer">Copyright 1987-2005 C. de Boor and The MathWorks, Inc.<br>            Published with MATLAB&reg; 7.1<br></p>      </div>      <!--##### SOURCE BEGIN #####%% Splines and real world data% This a demonstration of a spline function being fit to real-world data.% Copyright 1987-2005 C. de Boor and The MathWorks, Inc.% $Revision: 5.24.4.2 $%% Manual knot choice for spline interpolation% Here are some data which record a property of% titanium measured as a function of temperature. We'll use it to% illustrate some issues with spline interpolation.%% The plot of the data shows a rather sharp peak.hold off[xx,yy]=titanium;frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];plot(xx,yy,'x');axis(frame);%%% We pick a few data points from these somewhat% rough data, since we want to interpolate. Here is a% picture of the data, with the selected data points marked.hold onpick=[1 5 11 21 27 29 31 33 35 40 45 49];tau=xx(pick);y=yy(pick);plot(tau,y,'ro');%%% Since a spline of order k with n+k knots has n degrees of% freedom, and we have 12 data points, a fit with a fourth order% spline requires 12+4 = 16 knots. Moreover, this knot sequence% t must be such that the i-th data site lies in the support% of the i-th B-spline. We achieve this by using the data% sites as knots, but adding two simple knots at either end.%%   dl = tau(2)-tau(1); dr = tau(end)-tau(end-1);%   t = [tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];%% We use this knot sequence to construct an interpolating cubic spline:%%   sp = spapi(t, tau, y);dl=tau(2)-tau(1);dr=tau(end)-tau(end-1);t=[tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]];   % construct the knot sequencehold onaxis(frame+[-2*dl 2*dr 0 0])plot(t,repmat(frame(3)+.03,size(t)),'kx')hold offsp=spapi(t,tau,y);                         % This constructs the spline.%%% Now, for the plot.% Since we do not care about the part of the spline% outside the data interval, we restrict the plot to that interval:%%  hold on,  fnplt(sp,[tau(1) tau(end)], 'k'),  hold off%plot(xx,yy,'x',tau,y,'ro'), axis(frame), hold on% Now, for the plot:fnplt(sp,[tau(1) tau(end)], 'k')hold off%%% A closer look at the left part of the spline fit shows some% undulations.xxx = linspace(tau(1),tau(5),41);plot(xxx, fnval(sp, xxx), 'k', tau, y, 'ro');axis([tau(1) tau(5) 0.6 1.2]);%%% The unreasonable bump in the first interval stems from the fact% that our spline goes smoothly to zero at its first knot.%% Here is a picture of the entire spline, along with its knot sequence.axis([tau(1) tau(5) 0.6 1.2]);fnplt(sp,'k');hold on, plot(t,repmat(.1,size(t)),'kx'), hold off%%% Here are again the data points as well.hold onplot(tau, y, 'ro');hold off%%% Here is a simple way to enforce a more reasonable boundary% behavior. We add two more data points outside the given data% interval and choose as our data there the values of the straight% line through the first two data points.tt=[tau(1)-[4 3 2 1]*dl tau tau(end)+[1 2 3 4]*dr];xx=[tau(1)-[2 1]*dl tau tau(end)+[1 2]*dr];yy=[y(1)-[2 1]*(y(2)-y(1)) y y(end)+[1 2]*(y(end)-y(end-1))];sp2=spapi(tt,xx,yy); fnplt(sp2,'b',tau([1 end]))hold onplot(tau,y,'or', xx([1 2 end-1 end]),yy([1 2 end-1 end]),'ko');axis(frame+[-2*dl 2*dr 0 0]);hold off%%% Here is also the original spline fit, shown in black, to show the% reduction of the undulation in the first and last interval.hold onfnplt(sp,'k',tau([1 end]))hold off%%% Finally, here is a closer look at the first four data intervals% which shows more clearly the reduction of the undulation near% the left end.plot(xxx,fnval(sp2,xxx),'b',tau,y,'ro',xxx,fnval(sp,xxx),'k');axis([tau(1) tau(5) .6 2.2]);hold off%% Automatic knot choice for interpolation% If all this detail turns you off, let the spline toolbox choose the knots% for you, by using the spline interpolation command SPAPI in the form%        spapi( order, data_sites, data_values )autosp = spapi( 4, tau, y);fnplt(autosp,'g')hold on, plot(tau, y, 'or'), hold off%%% Here is the result of a much better knot choice (plotted above as red x's),% obtained by shifting the knot at 842 slightly to the right% and the knot at 985 slightly to the left.knots = fnbrk(autosp,'knots');hold on, plot(knots, repmat(.5,size(knots)),'xg')knots([7 12]) = [851, 971];plot(knots, repmat(.54,size(knots)),'xr')adjsp = spapi(knots, tau, y);fnplt(adjsp,'r',2), hold off%%% Else, simply try the standard cubic spline interpolant,% supplied by CSAPI (which amounts to a slightly different choice of knots):%% The next and last slide shows all five interpolants, for comparison.autocs = csapi(tau, y);fnplt(autocs,'c')hold on, plot(tau, y, 'or'), hold off%%% With such fast varying data, it is hard to get agreement among all reasonable% interpolants, even if each of them is a cubic spline.hold onfnplt(sp,'k',tau([1 end]))  % black: originalfnplt(sp2,'b',tau([1 end])) % blue:  with special end conditionsfnplt(autosp,'g')           % green: automatic knot choice by SPAPIfnplt(autocs,'c')           % cyan:  automatic knot choice by CSAPIfnplt(adjsp,'r',2)          % red:   knot choice by SPAPI slightly changedhold offdisplayEndOfDemoMessage(mfilename)##### SOURCE END #####-->   </body></html>
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -