多项式曲线拟合

1121865306 · 发表于 2016-4-25 13:34:46

多项式曲线拟合，求助。
有一组数据，想求得多项式拟合的表达式。
有没有编好的，或者IMSL里多项式拟合的方法。
自己写好难，也比较着急。
跪谢！！

pasuka · 发表于 2016-4-25 20:16:43

短平快的话，为啥不考虑Excel或者MATLAB呢？

1121865306 · 发表于 2016-4-25 21:42:42

后面还要做别的，所以打算统一都用fortran了，

1121865306 · 发表于 2016-4-25 21:42:58

已经找到例子了，谢谢平台

taipinghu · 发表于 2017-10-7 16:43:46

1121865306 发表于 2016-4-25 21:42
已经找到例子了，谢谢平台

楼主方便传阅一下例子吗？最近在搞这个问题。谢谢楼主！

kerb · 发表于 2017-11-20 15:58:10

本帖最后由 kerb 于 2018-7-22 14:18 编辑

这是我写的一个，希望对你有用

全域拉格朗日插值其实就是多项式插值，只是一般教科书的讲解形式，让很多人忽略了这一点。高阶多项式插值会引起容格现象，也就是说虽然多项式插值会经过已知点，但是会引起数值震荡。
最小二乘法适用于数据拟合，最常见的是那种线性拟合，其实也可以用已知的曲线的组合去拟合给定的数据。曲线可以有多种选择，比如光谱曲线往往用高斯曲线或者洛仑兹曲线的组合去拟合，当然也可以用多项式曲线去拟合，只是缺少了物理意义，但是这个时候，可以把这种拟合的曲线理解为插值曲线，给出数值区间任意点的函数值。

问题描述：已知一组数据点(XT,YT)，这些点的总数是M，我们希望用一个N阶的多项式去拟合这些数据，假设这个多项式的系数用数组U表示，记住一个N阶的多项式，有N+1个系数（包括常数项），求的多项式系数，多项式就确定了。

Y(1)=U(1)+U(2)*X(1)+U(3)*X(1)^2+...+U(N+1)*X(1)^N
Y(2)=U(1)+U(2)*X(2)+U(3)*X(2)^2+...+U(N+1)*X(2)^N
......................................................................................
Y(M)=U(1)+U(2)*X(M)+U(3)*X(M)^2+...+U(N+1)*X(M)^N
记
Y=(Y(1),Y(2),...,Y(M))^T，U=(U(1),U(2),...,U(N))
矩阵A的元素是A(i,j)=X(i)^j，1<=i<=M，0<=j<=N，A是MxN的
上面的方程组可以写成：
AU=Y
如果M<N，这个时候有无穷多解，但是无法解释每个解的物理含义；如果M=N，且A非奇异，这个时候可以解得一组数据，这个时候其实就是求多项式插值；当M>N，这是超定方程组，理论上只能求得最小二乘解，记A^T表示矩阵A的转置矩阵
(A^T*A)U=(A^T*Y)，解出U就求出多项式系数。

下面的程序中：
POINTS_NUM表示已知点的数目，即上面说的M
POLY_ORDER,表示多项式阶数，即上面说的N
POLY_COEF多项式系数，数组
XT，YT已知点，数组

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

001PROGRAM MAIN !主程序
002    IMPLICIT NONE
003    INTEGER,PARAMETER::POINTS_NUM=15 !这个是测试点的数目，你应该根据你的测试点数修改它
004    INTEGER,PARAMETER::POLY_ORDER=13 !如果是线性拟合：POLY_ORDER=1
005    INTEGER,PARAMETER::DP=KIND(1.D0)
006    REAL(KIND=DP)::A(POINTS_NUM,POLY_ORDER+1),POLY_COEF(POLY_ORDER+1) !如果是线性拟合，实际上就是POLY_COEF(2)：y = ax + b
007    REAL(KIND=DP)::XT(POINTS_NUM),YT(POINTS_NUM),DX,XN(POINTS_NUM),YN(POINTS_NUM)
008    INTEGER::INFO
009    INTEGER::I
010   !XT,YT为测试的数据
011    XT=(/0.1D0,0.2D0,0.3D0,0.4D0,0.5D0,0.6D0,0.7D0,0.8D0,0.9D0,0.92D0,0.94D0,0.96D0,0.98D0,0.99D0,0.995D0/)!插值节点
012    YT=(/0.86D0,0.78D0,0.71D0,0.66D0,0.59D0,0.52D0,0.45D0,0.36D0,0.23D0,0.19D0,0.16D0,0.12D0,0.07D0,0.04D0,0.02D0/)!插值节点函数值
013    !XN=(/0.15D0,0.985D0,0.995D0/)!待计算的插值点
014    !YN=(/0.86060402972635042,5.30099995722376274E-002,2.57654880479794013E-002/)!插值结果
015    DX=1.D0/15
016    CALL LSINTERPOLY(POINTS_NUM,XT,YT,POLY_ORDER,POLY_COEF,INFO)
017    WRITE(*,*)
018    WRITE(*,*)
019    IF(INFO/=-1)THEN
020        DO I=1,15
021            XN(I)=XT(I)
022            CALL POLY_VALUE(POLY_ORDER,POLY_COEF,XN(I),YN(I))
023            !WRITE(*,FMT="(4X,A,G0,A,G0)")"YN(",XN(I),")=",YN(I)
024            WRITE(*,FMT="(F6.4,2X,G0)")YN(I),YN(I)-YT(I)
025        END DO
026    ELSE
027        WRITE(*,*)"ERROR"
028    ENDIF
029    READ(*,*)
030END PROGRAM MAIN
031!---------------------------------------------------------------------------
032SUBROUTINE LSINTERPOLY(POINTS_NUM& !测试数据的个数
033                                          ,XT& !测试数据横坐标
034                                          ,YT& !测试数据纵坐标
035                                          ,POLY_ORDER& !多项式阶数
036                                          ,POLY_COEF& !多项式系数
037                                          ,INFO& !输出信息，INF=0，求解成功，INFO=-1求解失败
038                                          )!求多项式系数
039    IMPLICIT NONE
040    INTEGER,PARAMETER::DP=KIND(1.D0)
041    INTEGER,INTENT(IN)::POINTS_NUM,POLY_ORDER
042    REAL(KIND=DP),INTENT(IN)::XT(POINTS_NUM),YT(POINTS_NUM)
043    REAL(KIND=DP),INTENT(OUT)::POLY_COEF(POLY_ORDER+1)
044    INTEGER,INTENT(OUT)::INFO
045    LOGICAL::SOL_OK
046    REAL(KIND=DP)::A(POINTS_NUM,POLY_ORDER+1),ATA(POLY_ORDER+1,POLY_ORDER+1)
047    REAL(KIND=DP)::B(POINTS_NUM),C(POLY_ORDER+1,POINTS_NUM)
048    INTEGER::I,J
049    IF(POLY_ORDER+1>POINTS_NUM)THEN
050        INFO=-1
051        RETURN
052    END IF
053    INFO=0
054    A=0.D0
055    DO I=1,POINTS_NUM
056        DO J=0,POLY_ORDER
057            IF(J==0)THEN
058                A(I,J+1)=1.D0
059            ELSE
060                A(I,J+1)=XT(I)**J
061            ENDIF
062        END DO
063    END DO
064    B=YT
065    C=TRANSPOSE(A)
066    ATA=MATMUL(C,A)
067    POLY_COEF=MATMUL(C,B)
068    
069!    IF(POLY_ORDER+1==POINTS_NUM)THEN
070!        CALL SOLVE(A,POLY_ORDER+1,POLY_COEF,SOL_OK)
071!    ELSE
072!        CALL SOLVE(ATA,POLY_ORDER+1,POLY_COEF,SOL_OK)
073!    END IF
074    CALL SOLVE(ATA,POLY_ORDER+1,POLY_COEF,SOL_OK)
075 
076    IF(SOL_OK)THEN
077        INFO=0
078    ELSE
079        INFO=-1
080    END IF
081    RETURN
082END SUBROUTINE
083!---------------------------------------------------------------------------
084SUBROUTINE POLY_VALUE(POLY_ORDER,POLY_COEF,X,Y)!根据多项式计算插值点(X)的函数值Y
085    IMPLICIT NONE
086    INTEGER,PARAMETER::DP=KIND(1.D0)
087    INTEGER,INTENT(IN)::POLY_ORDER
088    REAL(KIND=DP),INTENT(IN)::POLY_COEF(POLY_ORDER+1),X
089    REAL(KIND=DP),INTENT(OUT)::Y
090    INTEGER::I
091    Y=0.D0
092    DO I=0,POLY_ORDER
093        Y=Y*X+POLY_COEF(POLY_ORDER+1-I)
094    END DO
095    RETURN
096END SUBROUTINE
097!---------------------------------------------------------------------------
098SUBROUTINE SOLVE(A,N,B,SOL_OK)!全选主元高斯消元回带法求解器
099    IMPLICIT NONE
100    INTEGER,PARAMETER::DP=KIND(1.D0)
101    INTEGER,INTENT(IN)::N
102    LOGICAL,INTENT(OUT)::SOL_OK
103    REAL(KIND=DP),INTENT(INOUT)::A(N,N)
104    REAL(KIND=DP),INTENT(INOUT)::B(N)
105    REAL(KIND=DP)::BIG,DUM,PIVINV,TEMP
106    INTEGER I,ICOL,IROW,J,K,L,LL,INDXC(N),INDXR(N),IPIV(N)
107    SOL_OK=.TRUE.
108    DO J=1,N
109        IPIV(J)=0
110    ENDDO
111    DO I=1,N
112        BIG=0.D0
113        DO J=1,N
114            IF(IPIV(J).NE.1)THEN
115                DO K=1,N
116                    IF (IPIV(K).EQ.0) THEN
117                        IF (ABS(A(J,K)).GE.BIG)THEN
118                            BIG=ABS(A(J,K))
119                            IROW=J
120                            ICOL=K
121                        ENDIF
122                    ELSEIF (IPIV(K).GT.1) THEN
123                        WRITE(*,*)'SINGULAR MATRIX IN GAUSSJ'
124                        SOL_OK=.FALSE.
125                    ENDIF
126                ENDDO
127            ENDIF
128        ENDDO
129        IPIV(ICOL)=IPIV(ICOL)+1
130        IF (IROW.NE.ICOL) THEN
131            DO L=1,N
132                DUM=A(IROW,L)
133                A(IROW,L)=A(ICOL,L)
134                A(ICOL,L)=DUM
135            ENDDO
136            DUM=B(IROW)
137            B(IROW)=B(ICOL)
138            B(ICOL)=DUM
139        ENDIF
140        INDXR(I)=IROW
141        INDXC(I)=ICOL
142        TEMP=ABS(A(ICOL,ICOL))
143        IF(A(ICOL,ICOL)==0.D0)THEN
144            WRITE(*,*) 'SINGULAR MATRIX IN GAUSSJ'
145            SOL_OK=.FALSE.
146            RETURN
147        ENDIF
148        PIVINV=1.D0/A(ICOL,ICOL)
149        A(ICOL,ICOL)=1.D0
150        DO L=1,N
151            A(ICOL,L)=A(ICOL,L)*PIVINV
152        ENDDO
153        B(ICOL)=B(ICOL)*PIVINV
154        DO LL=1,N
155            IF(LL.NE.ICOL)THEN
156                DUM=A(LL,ICOL)
157                A(LL,ICOL)=0.D0
158                DO L=1,N
159                    A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
160                ENDDO
161                B(LL)=B(LL)-B(ICOL)*DUM
162            ENDIF
163        ENDDO
164    ENDDO
165    DO L=N,1,-1
166        IF(INDXR(L).NE.INDXC(L))THEN
167            DO K=1,N
168                DUM=A(K,INDXR(L))
169                A(K,INDXR(L))=A(K,INDXC(L))
170                A(K,INDXC(L))=DUM
171            ENDDO
172        ENDIF
173    ENDDO
174    RETURN
175END SUBROUTINE SOLVE

渡箭 · 发表于 2018-5-30 19:42:50

很给力！

rovingboy · 发表于 2018-6-26 22:14:22

DO L=N,1,-1
      IF(INDXR(L).NE.INDXC(L))THEN
         DO K=1,N
         DUM=A(K,INDXR(L))
         A(K,INDXR(L))=A(K,INDXC(L))
         A(K,INDXC(L))=DUM
      ENDDO
      ENDIF
ENDDO

rovingboy · 发表于 2018-6-26 22:15:28

其他都看懂了，最后一段的上面这个代码不知道意义何在？这个时候的A不是除了对角线为1，其他都为0吗？而且这个A值也没有返回其他程序里面，为什么需要这一段呢？

luott_33 · 发表于 2019-11-14 22:18:23

kerb 发表于 2017-11-20 15:58
这是我写的一个，希望对你有用

全域拉格朗日插值其实就是多项式插值，只是一般教科书的讲解形式，让很多人 ...

067行少了一步矩阵求逆？

[插值拟合] 多项式曲线拟合

评分

水王勋章

元老勋章

热心勋章

规矩勋章

爱心勋章