三次样条插值就是分段的低次（3次）插值多项式

总结：当插值点较多的时候，全局拉格朗日方法会出现插值不稳定，差值结果误差很大的情况；后来采用了三次样条函数插值的方法，差值效果很好，解决了插值发散的问题。

发现拉格朗插值的稳定性不好，下面的数据插值后发散。期待雪球用三次样条验证下是否可行（我改了下fortran代码，老是不对）？
共53行数据，要求插值间距为0.1：

0 0.0448691
0.8 0
3.36 -0.141092
5.92 -0.270408
8.48 -0.379768
11.04 -0.462291
13.6 -0.511338
16.16 -0.524469
18.72 -0.509841
21.28 -0.47633
23.84 -0.432552
26.4 -0.387011
28.96 -0.346568
31.52 -0.311683
34.08 -0.281601
36.64 -0.255941
39.2 -0.234728
41.76 -0.218125
44.32 -0.204685
46.88 -0.192656
49.44 -0.180615
52 -0.16741
54.56 -0.152814
57.12 -0.137541
59.68 -0.122244
62.24 -0.107776
64.8 -0.0951893
67.36 -0.0854269
69.92 -0.0777986
72.48 -0.071234
75.04 -0.0648455
77.6 -0.0579076
80.16 -0.0502016
82.72 -0.0420509
85.28 -0.0337848
87.84 -0.0258565
90.4 -0.0188259
92.96 -0.0132322
95.52 -0.00877765
98.08 -0.00491217
100.64 -0.0011118
103.2 0.0032074
105.76 0.00809623
108.32 0.0127307
110.88 0.0166258
113.44 0.019518
116 0.0215117
118.56 0.0226838
121.12 0.0216095
123.68 0.0171785
126.24 0.00933951
128.8 0
129.6 -2.70E-03

fcode 发表于 2014-1-24 15:26
在一定已知情况的前提下，调整插值算法，可以得到更符合真实的结果。

所以要具体问题具体分析。 ...

嗯，具体问题具体分析

在一定已知情况的前提下，调整插值算法，可以得到更符合真实的结果。

所以要具体问题具体分析。

谢谢！这个很完美，下午看。

为什么不用三次样条呢？我觉得效果还不错。

[Fortran] 纯文本查看 复制代码

Program www_fcode_cn
  use SplineMod
  Implicit None
  Integer , parameter :: N = 10
  Integer , parameter :: M = 37
  Real :: x(N) , y(N) , xx(M) , yy(M)
  integer :: i
  Open( 12 , File = 'a.txt' )
  Open( 13 , File = 'b.txt' )
  Do i = 1 , N
    Read( 12 , * ) x(i) , y(i)
  End Do
  call SplineMod_Init( N, x , y )
  Do i = 1 , M
    xx(i) = (i-1) * 0.5
    yy(i) = SplineMod_Interp( N , x , y , xx(i) )
  End Do
  Do i = 1 , M
    Write( 13 , * ) xx(i) , yy(i)
  End Do
End Program www_fcode_cn

Module SplineMod
  Implicit None  
  Integer , private ::  i , N
  Real , private , Allocatable :: S2(:) , Dely(:) , S3(:)

Contains

  Subroutine SplineMod_Init( iN , x , y )
    Integer , Intent( IN )  :: iN
    Real , Intent( IN )      :: x( iN ) , y( iN )
    Real :: B( iN ) , C( iN ) , H( iN ) , H1( iN ) , Delsqy
    Integer jj , N1
    N  = iN
    N1 = N - 1
    Allocate( S2( N ) )
    Allocate( Dely( N ) )
    Allocate( S3( N ) )
    Do i = 1 , N1
      H( i ) = x( i+1 ) - x( i )
      if ( abs(H(i)) < 1./3600 ) H(i) = 1./3600
      Dely( i ) = ( y(i+1) - y(i) ) / H( i ) 
    End Do
    Do i = 2 , N1 
      H1( i ) = H( i-1 ) + H( i ) 
      B( i )  = 0.5 * H(i-1) / H1( i ) 
      Delsqy = ( Dely(i) - Dely(i-1) ) / H1( i ) 
      S2( i ) = 2.0 * Delsqy
      C( i ) = 3.0 * Delsqy
    End Do
    S2( 1 ) = 0.0 
    S2( N ) = 0.0 
    Do jj = 1 , 26 
      Do i = 2 , N1 
        S2(i) = (C(i)-B(i)*S2(i-1)-(0.5-B(i))*S2(i+1)-S2(i))*1.0717968+S2(i) 
      End Do
    End Do
    Do i = 1 , N1 
      S3( i ) = ( S2(i+1) - S2(i) ) / H( i ) 
    End Do
  End Subroutine SplineMod_Init

  Real Function SplineMod_Interp( iN , x , y , T )
    Integer , Intent( IN ) :: iN
    Real , Intent( IN ) :: T
    Real , Intent( IN )      :: x( iN ) , y( iN )
    Integer i
    Real                 :: ht1 , ht2 , Delsqs    
    i = 1 
    if( ( T - x(i) ) <= 0.0 ) goto 17 
    if( ( T - x(N) ) <  0.0 )  goto 57 
    goto 59 
 56 if( ( T - x(i) ) <  0.0 ) goto 60 
    if( ( T - x(i) ) == 0.0 ) goto 17 
 57 i = i + 1 
    GOTO 56 
 59 i = N 
 60 i = i - 1 
 17 HT1 = T - x(i)
    HT2 = T - x(i+1)
    Delsqs = ( 2.0 * S2(i) + S2(i+1) + HT1 * S3(i) ) / 6.0 
    SplineMod_Interp = y(i) + HT1 * Dely( i ) + HT1 * HT2 * Delsqs 
  End Function SplineMod_Interp

  Subroutine SplineMod_UnInit()
    DeAllocate( S2 )
    DeAllocate( Dely )
    DeAllocate( S3 )
  End Subroutine SplineMod_UnInit
  
End Module SplineMod