【原创】Hermite插值的实现及各段系数的提取

dynamic

在Matlab中，我们可以使用pwch()函数来实现

p=pwch(x,y,s)返回插值系数，我们可以使用ppval(p,x0)来计算某一点的插值结果

但是请记住pwch返回的不是每一段的直接多项式的系数，例如下面的表达错误的

>> x = 1:10;
>> y = sin( x );
>> s = cos( x ); % slopes
>> pp = pwch( x, y, s ) ；%interpolation
>>polyval(pp.coefs(1,:),x(1))%这个用法是错误的

而实际上pwch返回的pp.coefs的第i行，是下面多项式的系数

f=coefs(i,1)*(X-breaks(i))^3 +coefs(i,2)*(X-breaks(i))^2 +coefs(i,3)*(X-breaks(i)) + coefs(i,4)

其中，i取从1到length(x)-1

Matlab只提供了怎么计算X对应插值的函数ppval()，没有提供直接多项式系数提取函数。也就是说,你要计算某些点的对应插值时，只需ppval(pp,xx)，这里pp是pwch函数的输出pp形式，必须照写，xx可以是标量或矩阵

而你是想得到形如c(i,1)*(X)^3 +c(i,2)*(X)^2 +c(i,3)*X + co(i,4)的各段上的多项式系数，就要把最上面的多项式 f  展开后，提取系数，而对此matlab没有提供直接可用的函数

下面的代码可以得到你需要的系数：

1. %by dsj6700417
   
2. %all rights reserved by www.matlabsky.cn
   
3. x = 1:10;
   
4. y = sin( x );
   
5. s = cos( x ); % slopes
   
6. pp = pwch( x, y, s ) ;
   
7. [breaks,coefs,l,k,d] = unmkpp(pp);
   
8. syms X,
   
9. c=zeros(9,4);
   
10. for i=1length(x)-1)
    
11. f=coefs(i,1)*(X-breaks(i))^3 +coefs(i,2)*(X-breaks(i))^2 ...
    
12. +coefs(i,3)*(X-breaks(i)) + coefs(i,4);
    
13. c(i,:) = sym2poly(f);
    
14. end
    
15. c  %展开后的插值多项式的系数

复制代码

注：如果你要计算x=1时，插值多项式的值，用polyval(pp.coefs(1,:),x(1))是不正确的
按照上面的程序，可以用下面两个办法得到：

1. >> ppval(pp,1)
   
2. ans =
   
3.     0.8415
   
4. >> polyval(c(1,:),1)
   
5. ans =
   
6.     0.8415

复制代码

再注：下面的代码更直接的得到了你想要的系数：

1. %by dsj6700417
   
2. %all rights reserved by www.matlabsky.cn
   
3. x = 1:10;
   
4. y = sin( x );
   
5. z = cos( x );
   
6. pp=zeros(length(x)-1,4);
   
7. for i=1length(x)-1)
   
8.     b=[y(i) y(i+1) z(i) z(i+1)]';
   
9.     A1=[x(i)^3 x(i)^2 x(i) 1];
   
10.     A2=[x(i+1)^3 x(i+1)^2 x(i+1) 1];
    
11.     A3=[3*x(i)^2 2*x(i) 1 0];
    
12.     A4=[3*x(i+1)^2 2*x(i+1) 1 0];  
    
13.     A=[A1; A2; A3; A4];
    
14.     pp(i,:)=inv(A)*b;
    
15.    end
    
16. pp　　%展开后的插值多项式的系数

复制代码

freeok

看不懂啊，，基础太差了。。