Fortran Coder

查看: 12707|回复: 12

[插值拟合] 关于差值方法

[复制链接]

98

帖子

5

主题

3

精华

专家

F 币
426 元
贡献
275 点

管理勋章新人勋章

发表于 2014-1-23 20:05:41 | 显示全部楼层 |阅读模式
已知一组数据,用Lagrange多项式插值方法计算出来的数据不符合数据的真实趋势,如何解决?

1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-1-23 20:07:31 | 显示全部楼层
数据算例拿出来瞧瞧

135

帖子

15

主题

0

精华

版主

F 币
1159 元
贡献
637 点

爱心勋章管理勋章

发表于 2014-1-23 22:38:07 | 显示全部楼层
上数据,Hermite插值

98

帖子

5

主题

3

精华

专家

F 币
426 元
贡献
275 点

管理勋章新人勋章

 楼主| 发表于 2014-1-24 08:43:23 | 显示全部楼层
这是已知数据:
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
下图是插值前后的对比图,显然插值后数据趋势不正确(黄色点为原始数据)。

插值前后对比图

插值前后对比图



1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-1-24 10:46:27 | 显示全部楼层
为什么不用三次样条呢?我觉得效果还不错。

PhoXo1.jpg

[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

98

帖子

5

主题

3

精华

专家

F 币
426 元
贡献
275 点

管理勋章新人勋章

 楼主| 发表于 2014-1-24 12:10:21 | 显示全部楼层
谢谢!这个很完美,下午看。

98

帖子

5

主题

3

精华

专家

F 币
426 元
贡献
275 点

管理勋章新人勋章

 楼主| 发表于 2014-1-24 15:16:41 | 显示全部楼层
再次感谢臭石头!
下午考虑的情况:三次样条虽然能很好地光滑插值(如臭石头给出的例子),但是和影响线的真实情况有所区别(如弯矩影响线有尖点),所以最终采取数据分段的方法,依旧采用Lagrange差值计算,结果吻合良好:
捕获.PNG

1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2014-1-24 15:26:50 | 显示全部楼层
在一定已知情况的前提下,调整插值算法,可以得到更符合真实的结果。

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

136

帖子

3

主题

0

精华

版主

F 币
1964 元
贡献
1677 点

帅哥勋章管理勋章爱心勋章新人勋章热心勋章元老勋章

发表于 2014-2-17 17:26:41 | 显示全部楼层
fcode 发表于 2014-1-24 15:26
在一定已知情况的前提下,调整插值算法,可以得到更符合真实的结果。

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

嗯,具体问题具体分析

98

帖子

5

主题

3

精华

专家

F 币
426 元
贡献
275 点

管理勋章新人勋章

 楼主| 发表于 2014-3-12 09:13:47 | 显示全部楼层

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

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

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-29 22:27

Powered by Tencent X3.4

© 2013-2024 Tencent

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