关于差值方法
已知一组数据，用Lagrange多项式插值方法计算出来的数据不符合数据的真实趋势，如何解决？数据算例拿出来瞧瞧 上数据，Hermite插值 这是已知数据：
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
下图是插值前后的对比图，显然插值后数据趋势不正确（黄色点为原始数据）。
为什么不用三次样条呢？我觉得效果还不错。
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
谢谢！这个很完美，下午看。 再次感谢臭石头！
下午考虑的情况：三次样条虽然能很好地光滑插值（如臭石头给出的例子），但是和影响线的真实情况有所区别（如弯矩影响线有尖点），所以最终采取数据分段的方法，依旧采用Lagrange差值计算，结果吻合良好：
在一定已知情况的前提下，调整插值算法，可以得到更符合真实的结果。
所以要具体问题具体分析。 fcode 发表于 2014-1-24 15:26
在一定已知情况的前提下，调整插值算法，可以得到更符合真实的结果。
所以要具体问题具体分析。 ...
嗯，具体问题具体分析
发现拉格朗插值的稳定性不好，下面的数据插值后发散。期待雪球用三次样条验证下是否可行（我改了下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
页:
[1]
2