关于差值方法

lm_lxt

已知一组数据，用Lagrange多项式插值方法计算出来的数据不符合数据的真实趋势，如何解决？

fcode

数据算例拿出来瞧瞧

珊瑚虫

上数据，Hermite插值

lm_lxt

这是已知数据：
0          0
2      0.350617
4       0.938272
6      0
8      -0.316049
10      -0.217284
12      0
14      0.0987654
16      0.0790123
18      0
下面是插值后加密的数据：
0              0
            0.5       -1.14289
              1      -0.970142
            1.5      -0.319806
              2       0.350617
            2.5       0.832135
              3        1.06738
            3.5        1.08199
              4       0.938272
            4.5       0.705968
              5       0.445728
            5.5       0.201669
              6              0
            6.5      -0.148687
              7      -0.245285
            7.5      -0.297569
              8      -0.316049
            8.5      -0.310836
              9      -0.289764
            9.5      -0.257765
             10      -0.217284
           10.5      -0.169407
             11      -0.115269
           11.5     -0.0573139
             12              0
           12.5      0.0503486
             13      0.0866936
           13.5       0.103369
             14      0.0987654
           14.5      0.0779355
             15      0.0540707
           15.5      0.0473299
             16      0.0790123
           16.5       0.158492
             17        0.25971
           17.5       0.283339
             18              0
下图是插值前后的对比图，显然插值后数据趋势不正确（黄色点为原始数据）。

fcode

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

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

lm_lxt

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

lm_lxt

再次感谢臭石头！
下午考虑的情况：三次样条虽然能很好地光滑插值（如臭石头给出的例子），但是和影响线的真实情况有所区别（如弯矩影响线有尖点），所以最终采取数据分段的方法，依旧采用Lagrange差值计算，结果吻合良好：

fcode

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

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

aliouying

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

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

嗯，具体问题具体分析

lm_lxt

发现拉格朗插值的稳定性不好，下面的数据插值后发散。期待雪球用三次样条验证下是否可行（我改了下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