查看: 9233|回复: 3|关注: 0

[已答复] 关于三次样条插值的默认边界条件 spline.m文件

[复制链接]

新手

5 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
发表于 2016-1-6 09:12:13 | 显示全部楼层 |阅读模式
最近在看三次样条插值,有需求将Matlab中的spline.m用高级语言重写,写到填充矩阵时,发现边界条件部分看不明白

else % set up the sparse, tridiagonal, linear system b = ?*c for the slopes
   b=zeros(yd,n);
   b(:,2:n-1)=3*(dx(dd,2:n-1).*divdif(:,1:n-2)+dx(dd,1:n-2).*divdif(:,2:n-1));
   if isempty(endslopes)
      x31=x(3)-x(1);xn=x(n)-x(n-2);
      b(:,1)=((dx(1)+2*x31)*dx(2)*divdif(:,1)+dx(1)^2*divdif(:,2))/x31;
      b(:,n)=...
      (dx(n-1)^2*divdif(:,n-2)+(2*xn+dx(n-1))*dx(n-2)*divdif(:,n-1))/xn;
   else
      x31 = 0; xn = 0; b(:,[1 n]) = dx(dd,[2 n-2]).*endslopes;
   end
   dxt = dx(:);
   c = spdiags([ [x31;dxt(1:n-2);0] ...
        [dxt(2);2*(dxt(2:n-1)+dxt(1:n-2));dxt(n-2)] ...
        [0;dxt(2:n-1);xn] ],[-1 0 1],n,n);



搜了百度发现这么一段,

在 Matlab中数据点称之为断点。 如果三次样条插值没有边界条件,最常用的方法,就是采用非扭结(not-a-knot)条件。
这个条件强迫第1个和第2个三次多项式的三阶导数相等。 对最后一个和倒数第 2个三次多项式也做同样地处理。这段话跟上面的代码是对应的吗,如何对应的啊,我没看明白,还请弄懂的指教一下。万分感谢,
回复主题 已获打赏: 0 积分

举报

新手

5 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2016-1-6 09:14:49 | 显示全部楼层
我已经看明白spline.m中是采用待定一阶导数的方法推导的三对角矩阵来解一阶导数值
但是怎么做到第1个和第2个三次多项式的三阶导数相等呢,在离散序列中,三阶导数怎么求啊?
回复此楼 已获打赏: 0 积分

举报

新手

5 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2016-1-6 09:19:55 | 显示全部楼层
spline.m的完整代码
% Check that data are acceptable and, if not, try to adjust them appropriately
[x,y,sizey,endslopes] = chckxy(x,y);
n = length(x); yd = prod(sizey);

% Generate the cubic spline interpolant in ppform

dd = ones(yd,1); dx = diff(x); divdif = diff(y,[],2)./dx(dd,:);
if n==2
   if isempty(endslopes) % the interpolant is a straight line
      pp=mkpp(x,[divdif y(:,1)],sizey);
   else         % the interpolant is the cubic Hermite polynomial
      pp = pwch(x,y,endslopes,dx,divdif); pp.dim = sizey;
   end
elseif n==3&&isempty(endslopes) % the interpolant is a parabola
   y(:,2:3)=divdif;
   y(:,3)=diff(divdif')'/(x(3)-x(1));
   y(:,2)=y(:,2)-y(:,3)*dx(1);
   pp = mkpp(x([1,3]),y(:,[3 2 1]),sizey);
else % set up the sparse, tridiagonal, linear system b = ?*c for the slopes
   b=zeros(yd,n);
   b(:,2:n-1)=3*(dx(dd,2:n-1).*divdif(:,1:n-2)+dx(dd,1:n-2).*divdif(:,2:n-1));
   if isempty(endslopes)
      x31=x(3)-x(1);xn=x(n)-x(n-2);
      b(:,1)=((dx(1)+2*x31)*dx(2)*divdif(:,1)+dx(1)^2*divdif(:,2))/x31;
      b(:,n)=...
      (dx(n-1)^2*divdif(:,n-2)+(2*xn+dx(n-1))*dx(n-2)*divdif(:,n-1))/xn;
   else
      x31 = 0; xn = 0; b(:,[1 n]) = dx(dd,[2 n-2]).*endslopes;
   end
   dxt = dx(:);
   c = spdiags([ [x31;dxt(1:n-2);0] ...
        [dxt(2);2*(dxt(2:n-1)+dxt(1:n-2));dxt(n-2)] ...
        [0;dxt(2:n-1);xn] ],[-1 0 1],n,n);

   % sparse linear equation solution for the slopes
   mmdflag = spparms('autommd');
   spparms('autommd',0);
   s=b/c;
   spparms('autommd',mmdflag);

   % construct piecewise cubic Hermite interpolant
   % to values and computed slopes
   pp = pwch(x,y,s,dx,divdif); pp.dim = sizey;

end

if nargin==2, output = pp; else output = ppval(pp,xx); end
回复此楼 已获打赏: 0 积分

举报

新手

5 麦片

财富积分


050


1

主题

4

帖子

0

最佳答案
发表于 2020-7-29 09:00:53 | 显示全部楼层
你好,就是我在在使用pp1 = csape(EGG_GPS,Q_1,'not-a-knot');时候遇到了问题,就是EGG_GPS和Q_1是对应的时间序列数据,为什么计算出的值不对,而且还出现了警告?
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表