注冊 登錄
微網社區 返回首頁

kerbcurb的個人空間 http://www.qgbobn.tw/?18294 [收藏] [復制] [分享] [RSS]

日志

多項式插值與最小二乘插值

已有 811 次閱讀2017-11-17 09:05 |系統分類:人文紀事

全域拉格朗日插值其實就是多項式插值,只是一般教科書的講解形式,讓很多人忽略了這一點。高階多項式插值會引起容格現象,也就是說雖然多項式插值會經過已知點,但是會引起數值震蕩。
最小二乘法適用于數據擬合,最常見的是那種線性擬合,其實也可以用已知的曲線的組合去擬合給定的數據。曲線可以有多種選擇,比如光譜曲線往往用高斯曲線或者洛侖茲曲線的組合去擬合,當然也可以用多項式曲線去擬合,只是缺少了物理意義,但是這個時候,可以把這種擬合的曲線理解為插值曲線,給出數值區間任意點的函數值。
 
問題描述:已知一組數據點(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表示已知點的數目
POLY_ORDER,表示多項式階數
POLY_COEF多項式系數數組
XT,YT已知點,數組
 
 
PROGRAM MAIN !主程序
    IMPLICIT NONE
    INTEGER,PARAMETER::POINTS_NUM=15
    INTEGER,PARAMETER::POLY_ORDER=13
    INTEGER,PARAMETER::DP=KIND(1.D0)
    REAL(KIND=DP)::A(POINTS_NUM,POLY_ORDER+1),POLY_COEF(POLY_ORDER+1)
    REAL(KIND=DP)::XT(POINTS_NUM),YT(POINTS_NUM),DX,XN(15),YN(15)
    INTEGER::INFO
    INTEGER::I
    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/)!插值節點
    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/)!插值節點函數值
    !XN=(/0.15D0,0.985D0,0.995D0/)!待計算的點
    !YN=(/0.86060402972635042,5.30099995722376274E-002,2.57654880479794013E-002/)
    DX=1.D0/15
    CALL LSINTERPOLY(POINTS_NUM,XT,YT,POLY_ORDER,POLY_COEF,INFO)
    WRITE(*,*)
    WRITE(*,*)
    IF(INFO/=-1)THEN
        DO I=1,15
            XN(I)=XT(I)
            CALL POLY_VALUE(POLY_ORDER,POLY_COEF,XN(I),YN(I))
            !WRITE(*,FMT="(4X,A,G0,A,G0)")"YN(",XN(I),")=",YN(I)
            WRITE(*,FMT="(F6.4,2X,G0)")YN(I),YN(I)-YT(I)
        END DO
    ELSE
        WRITE(*,*)"ERROR"
    ENDIF
    READ(*,*)
END PROGRAM MAIN
!---------------------------------------------------------------------------
SUBROUTINE LSINTERPOLY(POINTS_NUM,XT,YT,POLY_ORDER,POLY_COEF,INFO)!求多項式系數
    IMPLICIT NONE
    INTEGER,PARAMETER::DP=KIND(1.D0)
    INTEGER,INTENT(IN)::POINTS_NUM,POLY_ORDER
    REAL(KIND=DP),INTENT(IN)::XT(POINTS_NUM),YT(POINTS_NUM)
    REAL(KIND=DP),INTENT(OUT)::POLY_COEF(POLY_ORDER+1)
    INTEGER,INTENT(OUT)::INFO
    LOGICAL::SOL_OK
    REAL(KIND=DP)::A(POINTS_NUM,POLY_ORDER+1),ATA(POLY_ORDER+1,POLY_ORDER+1)
    REAL(KIND=DP)::B(POINTS_NUM),C(POLY_ORDER+1,POINTS_NUM)
    INTEGER::I,J
    IF(POLY_ORDER+1>POINTS_NUM)THEN
        INFO=-1
        RETURN
    END IF
    INFO=0
    A=0.D0
    DO I=1,POINTS_NUM
        DO J=0,POLY_ORDER
            A(I,J+1)=XT(I)**J
        END DO
    END DO
    B=YT
    C=TRANSPOSE(A)
    ATA=MATMUL(C,A)
    POLY_COEF=MATMUL(C,B)
   
!    IF(POLY_ORDER+1==POINTS_NUM)THEN
!        CALL SOLVE(A,POLY_ORDER+1,POLY_COEF,SOL_OK)
!    ELSE
!        CALL SOLVE(ATA,POLY_ORDER+1,POLY_COEF,SOL_OK)
!    END IF
    CALL SOLVE(ATA,POLY_ORDER+1,POLY_COEF,SOL_OK)

    IF(SOL_OK)THEN
        INFO=0
    ELSE
        INFO=-1
    END IF
    RETURN
END SUBROUTINE
!---------------------------------------------------------------------------
SUBROUTINE POLY_VALUE(POLY_ORDER,POLY_COEF,X,Y)!根據多項式計算插值點(X)的函數值Y
    IMPLICIT NONE
    INTEGER,PARAMETER::DP=KIND(1.D0)
    INTEGER,INTENT(IN)::POLY_ORDER
    REAL(KIND=DP),INTENT(IN)::POLY_COEF(POLY_ORDER+1),X
    REAL(KIND=DP),INTENT(OUT)::Y
    INTEGER::I
    Y=0.D0
    DO I=0,POLY_ORDER
        Y=Y*X+POLY_COEF(POLY_ORDER+1-I)
    END DO
    RETURN
END SUBROUTINE
!---------------------------------------------------------------------------
SUBROUTINE SOLVE(A,N,B,SOL_OK)!全選主元高斯消元回帶法求解器
    IMPLICIT NONE
    INTEGER,PARAMETER::DP=KIND(1.D0)
    INTEGER,INTENT(IN)::N
    LOGICAL,INTENT(OUT)::SOL_OK
    REAL(KIND=DP),INTENT(INOUT)::A(N,N)
    REAL(KIND=DP),INTENT(INOUT)::B(N)
    REAL(KIND=DP)::BIG,DUM,PIVINV,TEMP
    INTEGER I,ICOL,IROW,J,K,L,LL,INDXC(N),INDXR(N),IPIV(N)
    SOL_OK=.TRUE.
    DO J=1,N
        IPIV(J)=0
    ENDDO
    DO I=1,N
        BIG=0.D0
        DO J=1,N
            IF(IPIV(J).NE.1)THEN
                DO K=1,N
                    IF (IPIV(K).EQ.0) THEN
                        IF (ABS(A(J,K)).GE.BIG)THEN
                            BIG=ABS(A(J,K))
                            IROW=J
                            ICOL=K
                        ENDIF
                    ELSEIF (IPIV(K).GT.1) THEN
                        WRITE(*,*)'SINGULAR MATRIX IN GAUSSJ'
                        SOL_OK=.FALSE.
                    ENDIF
                ENDDO
            ENDIF
        ENDDO
        IPIV(ICOL)=IPIV(ICOL)+1
        IF (IROW.NE.ICOL) THEN
            DO L=1,N
                DUM=A(IROW,L)
                A(IROW,L)=A(ICOL,L)
                A(ICOL,L)=DUM
            ENDDO
            DUM=B(IROW)
            B(IROW)=B(ICOL)
            B(ICOL)=DUM
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        TEMP=ABS(A(ICOL,ICOL))
        IF(A(ICOL,ICOL)==0.D0)THEN
            WRITE(*,*) 'SINGULAR MATRIX IN GAUSSJ'
            SOL_OK=.FALSE.
            RETURN
        ENDIF
        PIVINV=1.D0/A(ICOL,ICOL)
        A(ICOL,ICOL)=1.D0
        DO L=1,N
            A(ICOL,L)=A(ICOL,L)*PIVINV
        ENDDO
        B(ICOL)=B(ICOL)*PIVINV
        DO LL=1,N
            IF(LL.NE.ICOL)THEN
                DUM=A(LL,ICOL)
                A(LL,ICOL)=0.D0
                DO L=1,N
                    A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
                ENDDO
                B(LL)=B(LL)-B(ICOL)*DUM
            ENDIF
        ENDDO
    ENDDO
    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
    RETURN
END SUBROUTINE SOLVE

路過

雞蛋

鮮花

握手

雷人

評論 (0 個評論)

facelist doodle 涂鴉板

您需要登錄后才可以評論 登錄 | 注冊

客服中心 搜索

QQ|友情鏈接||手機無線|微網社區 ( 蜀ICP備09016811號 )

GMT+8, 2020-2-29 02:53 , Processed in 0.062500 second(s), 18 queries .

Powered by Discuz! X3.4

Copyright © 2008-2020 Design: Comiis.Com

返回頂部
新浪排球