计算数学
数据拟合
多项式拟合星历
本节与下一节内容介绍用多项式和切比雪夫多项式拟合星历问题。实际计算中，由于数据较大，一般通过 逐行读取文件的方法。本节代码可以编译成两个程序，一个是拟合程序，生成拟合系数。一个是通过系数计算星历程序。 下一节内容用法同本节。拟合的Makefile文件如下。
objects = main.o add_neq.o basefunc.o \
chol_eq.o chol_rf.o covariance.o inv_dtri.o \
polyfit.o polyval.o ephfit.o ephfit_shell.f90
cc = ifort
co = -g -c
main.app:$(objects)
$(cc) -g -o ../../main.app $(objects)
main.o: main.f90
$(cc) $(co) main.f90
add_neq.o : add_neq.f90
$(cc) $(co) add_neq.f90
basefunc.o : basefunc.f90
$(cc) $(co) basefunc.f90
chol_eq.o: chol_eq.f90
$(cc) $(co) chol_eq.f90
chol_rf.o: chol_rf.f90
$(cc) $(co) chol_rf.f90
covariance.o: covariance.f90
$(cc) $(co) covariance.f90
inv_dtri.o: inv_dtri.f90
$(cc) $(co) inv_dtri.f90
ephfit.o: ephfit.f90
$(cc) $(co) ephfit.f90
polyfit.o: polyfit.f90
$(cc) $(co) polyfit.f90
polyval.o: polyval.f90
$(cc) $(co) polyval.f90
ephfit_shell.o: ephfit_shell.f90
$(cc) $(co) ephfit_shell.f90
objects = main.obj add_neq.obj basefunc.obj \
chol_eq.obj chol_rf.obj covariance.obj inv_dtri.obj \
polyfit.obj polyval.obj ephfit.obj ephfit_shell.obj
cc = ifort
co = -c
main.exe:$(objects)
$(cc) -o main.exe $(objects)
$(obj): %.obj: %.f90
$(cc) $(co) $< $@
subroutine add_neq(ATA,ATb,oc,H,N)
!---------------------------------subroutine comment
! Purpose : 2018-09-05 13:33 (Created)
! add a new observation to normal equations
! (ATA+ H'H) x = ATb + H'oc
! Changes :
! *.
!-----------------------------------------------------
! Input Parameters :
! ATA ---- initial normal matrix
! Atb ---- initial right vector
! oc ----- o-c for liner equation
! H ------ obs vector
! N ------- the dimention of the matrix
! Output Parameters :
! ATA ------ updated normal matrix
! ATb ------- updated right vector
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2018)
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8:: ATA(N,N),ATb(N) ,H(N)
do i = 1,N
if (H(i)==0d0) cycle
do j =1 ,N
ATA(i,j) = ATA(i,j)+ H(i)*H(j)
end do
ATb(i)= ATb(i)+H(i)*oc
end do
end subroutine add_neq
subroutine basefunc(bv,x,m)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date :2010.07.09
!-----------------------------------------------------
! Purpose : 任意阶多项式基函数
!
! Post Script :
! 1.
! 2.
! 3.
!-----------------------------------------------------
! Input parameters :
! 1. x 标量
! 2. m 多项式阶数
! Output parameters :
! 1.
! 2.
! Common parameters :
! 1.
! 2.
!----------------------------------------------------
implicit real*8(a-z)
integer::M
real*8::bv(m)
integer::i
bv(1)=1d0
do i=2,m
bv(i)=bv(i-1)*x
end do
end subroutine basefunc
subroutine chol_eq(A,b,x,P,N)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : song.yz
! Date : 2019.04.15
!-----------------------------------------------------
! Purpose : 对称正定方程解算（如法方程）
!
! Post Script :
! 1. 修正的不开平方Cholesky分解法
! 2.
! 3.
!-----------------------------------------------------
! Input parameters :
! 1. A(N,N)对称正定系数矩阵
! 2. b(N)右向量
! 3. N--方程维数
! 4.
! Output parameters :
! 1. x-计算结果
! 2. P ----协方差 （由于A已经为对称正定矩阵<法矩阵>，因此P为A的逆）
! 3.
! 4.
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit real*8(a-z)
integer::N
real*8::A(N,N),b(n),x(n),P(N,N)
!---------------------------------subroutine variable
real*8::L(N,N),d(n)
integer::i,k
real*8::y(N)
call chol_rf(A,L,d,P,n)
!以上已经得到Cholesky分解
y(1)=b(1)
do i=2,n
tmp1=0d0
do k=1,i-1
tmp1=tmp1+L(i,k)*y(k)
end do
y(i)=b(i)-tmp1
end do
x(n)=y(n)/d(n)
do i=n-1,1,-1
tmp1=0d0
do k=i+1,n
tmp1=tmp1+L(k,i)*x(k)
end do
x(i)=y(i)/d(i)-tmp1
end do
end subroutine chol_eq
subroutine chol_rf(A,L,d,P,N)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date : 2019.04.15
!-----------------------------------------------------
! Purpose :
!
! Post Script :
! 1. 不用开平方的Cholesky分解
! 2.
! 3.
!-----------------------------------------------------
! Input parameters :
! 1. A(N,N)--输入矩阵
! 2. N----矩阵维书
! 3.
! 4.
! Output parameters :
! 1. L矩阵
! 2. d---对角矩阵（用向量存储）
! 3.
! 4. P ---协方差矩阵（由于A已经为对称正定矩阵<法矩阵>，因此P为A的逆）
!-----------------------------------------------------
! Copyrigt :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit real*8(a-z)
integer::N
real*8::A(n,n),L(n,n),d(n),P(N,N)
!---------------------------------subroutine variable
integer::i,j
real*8:: g(n,n)
!设置初值
L=0d0
d(1)=a(1,1)
do i=2,n
do j=1,i-1
!此层循环算g(i,j)
tmp1=0d0
do k=1,j-1
tmp1=tmp1+g(i,k)*L(j,k)
end do
g(i,j)=a(i,j)-tmp1
end do
do j=1,i-1
!此层循环算L(i,j)
L(i,j)=g(i,j)/d(j)
end do
!以下算d(i)
tmp1=0d0
do k=1,i-1
tmp1=tmp1+g(i,k)*L(i,k)
end do
d(i)=a(i,i)-tmp1
end do
!设置对角线元素
do i=1,N
L(i,i)=1d0
end do
call covariance(L,D,P,N)
end subroutine chol_rf
subroutine covariance(L,D,P,N)
!-----------------------------------------------------
! Purpose : 2019-04-15 13:08 (Created)
! 根据LDL'分解计算协方差矩阵
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! L -----
! D ---- 对角线元素
! N ---- 矩阵维数
! Output Parameters :
!
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N
real*8 :: L(N,N),D(N),P(N,N)
!-----------------------------------------------------
integer :: i
real*8 :: invL(N,N),invLT(N,N)
!-----------------------------------------------------
call inv_dtri(L,invL,N)
invLT = transpose(invL)
do i =1,N
invLT(:,i)=invLT(:,i)/D(i)
end do
P = matmul(invLT,invL)
end subroutine covariance
subroutine ephfit(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:50 (Created)
!
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! M ----- dimention of observation data
! N ----- order of polynomial
!
! Output Parameters :
! px ---- coefficients of x
! py ---- coefficients of y
! pz ---- coefficients of z
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer ::i,ioerr,N, M ,iter
real*8 :: ATAx(N+1,N+1),ATbx(N+1),px(N+1),covx(N+1,N+1),Hx(N+1)
real*8 :: ATAy(N+1,N+1),ATby(N+1),py(N+1),covy(N+1,N+1),Hy(N+1)
real*8 :: ATAz(N+1,N+1),ATbz(N+1),pz(N+1),covz(N+1,N+1),Hz(N+1)
integer :: year,mon,day,hour,minu ,day0,hour0,minu0
real*8 :: sec,sec0, x, y ,z ,t, tspan
open(unit=11,file="xyz_ECF")
open(unit=17,file="tmp")
ATAx = 0D0
ATbx = 0d0
ATAy = 0D0
ATby = 0d0
ATAz = 0D0
ATbz = 0d0
iter = 0
do
read(11,*,iostat=ioerr) year,mon,day0,hour0,minu0,sec0 , x, y, z
if(ioerr/=0) exit
backspace(11)
iter = iter+1
write(12, '(A3,i5)') 'coe',iter
write(17, '(A3,i5)') 'coe',iter
do
read(11,*,iostat=ioerr) year,mon,day,hour,minu,sec , x, y, z
if(ioerr/=0) exit
!t = day - day0 + hour/24d0 + minu/1440d0 + sec/1440d0/3600d0
!t = (day - day0) * 24 - hour0 + hour + minu/60d0 + sec/3600d0
t= (day-day0)*24*60 + (hour-hour0)*60 + (minu - minu0) + (sec - sec0)/60d0
if (t>=tspan) then
backspace(11)
exit
end if
write(17,'(F15.8,I7,4I4,f10.3,3f20.3)') t, year,mon,day,hour,minu,sec, x, y, z
x = x/1d3
y = y/1d3
z = z/1d3
call basefunc(Hx,t,N+1)
call add_NEQ(ATAx,ATbx,x ,Hx,N+1)
call basefunc(Hy,t,N+1)
call add_NEQ(ATAy,ATby,y ,Hy,N+1)
call basefunc(Hz,t,N+1)
call add_NEQ(ATAz,ATbz,z ,Hz,N+1)
end do
call chol_eq(ATAx,ATbx,Px,covx,N+1)
call chol_eq(ATAy,ATby,Py,covy,N+1)
call chol_eq(ATAz,ATbz,Pz,covz,N+1)
ATAx = 0D0
ATbx = 0d0
ATAy = 0D0
ATby = 0d0
ATAz = 0D0
ATbz = 0d0
!write(12, '(A3,i5, I7,4I4,f10.3)') 'coe' ,iter , year,mon,day,hour,minu,sec
write(12,'((F25.12,/))') (px(i),i=1,N+1)
write(12,'((F25.12,/))') (py(i),i=1,N+1)
write(12,'((F25.12,/))') (pz(i),i=1,N+1)
end do
end subroutine ephfit
subroutine ephfit_shell(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-06-05 17:42 (Created)
!
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
!
!
!
! Output Parameters :
!
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N,i
real*8 :: px(N+1),py(N+1),pz(N+1)
real*8 :: tspan
!-----------------------------------------------------
open(unit=12,file='coe.txt')
write(12,'(A)') 'order of polynomial and end of time'
write(12,'(I5,F18.10)') N, tspan
call ephfit(N,tspan)
end subroutine ephfit_shell
subroutine inv_dtri(L,invL,N)
!-----------------------------------------------------
! Purpose : 2019-04-15 10:45 (Created)
! 计算下三角矩阵的逆矩阵，回带过程,无开平方运算
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! L ----- 下三角矩阵
! N ----- 矩阵维数
!
! Output Parameters :
! invL ------ 下三角的逆矩阵
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N
real*8 :: L(N,N),invL(N,N)
!----------------------------------------------------
integer :: i , j ,k
real*8 :: xsum
!-----------------------------------------------------
invL = 0d0
invL(1,1) = 1/L(1,1)
do i = 2, N
invL(i,i)= 1/L(i,i)
do j = 1,i-1
xsum = 0d0
do k = 1,i
xsum = xsum + L(i,k)*invL(k,j)
end do
invL(i,j)= - xsum / L(i,i)
end do
end do
end subroutine inv_dtri
program main
!-------------------------------------------- comment
! Purpose :
! Author : Song Yezhi
!
! Versions and Changes :
! V1.0 ----2019-05-08 09:16:37
!
!-----------------------------------------------------
! Notes :
! the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y.
! The coefficients in p are in powers
!
!
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit none
integer :: N, narg
character*8 :: strN,strSpan
real*8 :: tspan !unit : minute
narg = iargc()
!get the number of args input
if (narg<2) then
write(*,*) 'usage:'
write(*,*) 'main.app N tspan'
write(*,*) ' N is the order of the polynomial '
write(*,*) ' tspan is time span for data fit '
else
call getarg(1,strN)
call getarg(2,strSpan)
read(strN,*) N
read(strSpan,*) tspan
call ephfit_shell(N,tspan)
end if
end program main
subroutine polyfit(x,y,M,p,N)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:40 (Created)
!
! The coefficients in p (px =p(1) + p(2)x + p(3)x*x + .... )
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! x ---- 自变量数组
! y ---- 应变量数组
! M ----- 数据维数
! N ----- 多项式阶数
! Output Parameters :
! p --参数 P维数为多项式阶数+1
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N,M
real*8 :: x(M),y(M),p(N+1)
integer :: i
real*8 :: H(N+1),ATA(N+1,N+1),ATb(N+1),PX(N+1,N+1)
!-----------------------------------------------------
ATA =0D0
ATb =0d0
do i=1,M
call basefunc(H,x(i),N+1)
call add_NEQ(ATA,ATb,y(i),H,N+1)
end do
call chol_eq(ATA,ATb,Px,p,N+1)
end subroutine polyfitsubroutine polyval(p,N,x,fx)
!-----------------------------------------------------
! Purpose : 2019-05-08 10:11 (Created)
! polyval(p,x) returns the value of a polynomial of degree n evaluated at x.
! The input argument p is a vector of length n+1 whose elements are the coefficients in
! powers of the polynomial to be evaluated.
! fx = p(1) + p(2)x + p(3)x*x + ....
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! p 多项式系数
! N 多项式阶数
! x 自变量
! Output Parameters :
! fx 应变量
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N ,i
real*8 :: p(N+1),x ,fx ,H(N+1)
!-----------------------------------------------------
call basefunc(H,x,N+1)
fx =p(1)
do i=2,N+1
fx =fx + p(i) *H(i)
end do
end subroutine polyval
objects = main_eval.o basefunc.o polyval.o eval.o \
polyvel.o
cc = ifort
co = -g -c
main_eval.app:$(objects)
$(cc) -g -o ../../main_eval.app $(objects)
main_eval.o: main_eval.f90
$(cc) $(co) main_eval.f90
basefunc.o : basefunc.f90
$(cc) $(co) basefunc.f90
eval.o: eval.f90
$(cc) $(co) eval.f90
polyval.o: polyval.f90
$(cc) $(co) polyval.f90
polyvel.o: polyvel.f90
$(cc) $(co) polyvel.f90
subroutine basefunc(bv,x,m)
!---------------------------------subroutine comment
! Version : V1.0
! Coded by : syz
! Date :2010.07.09
!-----------------------------------------------------
! Purpose : 任意阶多项式基函数
!
! Post Script :
! 1.
! 2.
! 3.
!-----------------------------------------------------
! Input parameters :
! 1. x 标量
! 2. m 多项式阶数
! Output parameters :
! 1.
! 2.
! Common parameters :
! 1.
! 2.
!----------------------------------------------------
implicit real*8(a-z)
integer::M
real*8::bv(m)
integer::i
bv(1)=1d0
do i=2,m
bv(i)=bv(i-1)*x
end do
end subroutine basefunc
subroutine eval(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:50 (Created)
!
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! M ----- dimention of observation data
! N ----- order of polynomial
!
! Output Parameters :
! px ---- coefficients of x
! py ---- coefficients of y
! pz ---- coefficients of z
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: i,ioerr,N, kkk ,iter ,ibreak
real*8 :: px(N+1),py(N+1),pz(N+1)
integer :: year,mon,day,hour,minu ,day0,hour0,minu0
real*8 :: sec, sec0,fx, fy ,fz ,t, tspan
character*50 :: line
real*8 :: x,y,z,vx ,vy, vz
open(unit=12,file="xyz_ECF")
open(unit=13,file="xyz_com")
open(unit=14,file="xyz_com_tmp")
iter = 0
ibreak = 0
do
do
read(11,iostat=ioerr,fmt='(A)') line
if( ioerr/=0) then
ibreak = 1
! 设置ibreak 主要为跳出外层循环， 也可以直接采用goto 跳出外层循环
exit
end if
ibreak = 0
if (line(1:3)=='coe') exit
end do
if (ibreak==1) exit
do kkk = 1,N+1
read(11,*) px(kkk)
end do
do kkk = 1,N+1
read(11,*) py(kkk)
end do
do kkk = 1,N+1
read(11,*) pz(kkk)
end do
iter = iter+1
write(14,'(A3,I5)') 'coe',iter
write(14,'((F25.12,/))') (px(i),i=1,N+1)
write(14,'((F25.12,/))') (py(i),i=1,N+1)
write(14,'((F25.12,/))') (pz(i),i=1,N+1)
read(12,*,iostat=ioerr) year,mon,day0,hour0,minu0,sec0 , x, y, z
if(ioerr/=0) exit
backspace(12)
do
read(12,*,iostat=ioerr) year,mon,day,hour,minu,sec ,x,y,z
if(ioerr/=0) then
ibreak = 1
exit
end if
t= (day-day0)*24*60 + (hour-hour0)*60 + (minu - minu0) + (sec - sec0)/60d0
if (t>=tspan) then
backspace(12)
exit
end if
call polyval(px,N,t,fx)
call polyval(py,N,t,fy)
call polyval(pz,N,t,fz)
call polyvel(px,N,t,vx)
call polyvel(px,N,t,vy)
call polyvel(pz,N,t,vz)
fx=fx*1d3
fy=fy*1d3
fz=fz*1d3
!-------------------------
vx = vx*1d3 /60d0
vy = vy*1d3 /60d0
vz = vz*1d3 /60d0
!it is very important **************
!-------------------------
write(13,'(I5,4i4,f10.5,3f20.3,3f20.7,2f15.6,f15.3)') year,mon,day,hour,minu,sec, fx,fy,fz,vx ,vy ,vz ,0.0 ,0.0 ,0.0
write(14,'(F15.8,I5,4i4,f10.5,3f20.3,3f20.7,2f15.6,f15.3)') t, year,mon,day,hour,minu,sec, fx,fy,fz,vx ,vy ,vz ,0.0 ,0.0 ,0.0
end do
end do
close(12)
close(13)
end subroutine eval
program main_eval
!-------------------------------------------- comment
! Purpose :
! Author : Song Yezhi
! Song Yezhi 中科院上海天文台
! 天文地球动力学研究中心
! Versions and Changes :
! V1.0 ----2019-05-08 09:16:37
!
!-----------------------------------------------------
! Notes :
! the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y.
! The coefficients in p are in powers
!
!
!-----------------------------------------------------
! Copyrigt (C) :
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences
!----------------------------------------------------
implicit none
integer :: N
real*8 :: tspan
open(unit=11,file='coe.txt')
read(11,*)
read(11,*)N ,tspan
call eval(N,tspan)
close(11)
end program main_eval
subroutine polyval(p,N,x,fx)
!-----------------------------------------------------
! Purpose : 2019-05-08 10:11 (Created)
! polyval(p,x) returns the value of a polynomial of degree n evaluated at x.
! The input argument p is a vector of length n+1 whose elements are the coefficients in
! powers of the polynomial to be evaluated.
! fx = p(1) + p(2)x + p(3)x*x + ....
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
! p 多项式系数
! N 多项式阶数
! x 自变量
! Output Parameters :
! fx 应变量
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N ,i
real*8 :: p(N+1),x ,fx ,H(N+1)
!-----------------------------------------------------
call basefunc(H,x,N+1)
fx =p(1)
do i=2,N+1
fx =fx + p(i) *H(i)
end do
end subroutine polyval
subroutine polyvel(p,N,t,vel)
!-----------------------------------------------------
! Purpose : 2019-06-05 19:24 (Created)
!
!
! Changes :
!
!-----------------------------------------------------
! Input Parameters :
!
!
!
! Output Parameters :
!
!
!-----------------------------------------------------
! Author : Song Yezhi
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS
! (All rights reserved, 2019)
!----------------------------------------------------
implicit none
integer :: N ,i
real*8 :: p(N+1)
real*8 :: vel ,t
real*8 :: H(N+1)
!-----------------------------------------------------
H(1) = 0d0
H(2) = 1d0
H(3) = t
do i = 4,N+1
H(i) = H(i-1) *t
end do
!这里H仅仅包含多项式，不包含前面导数系数，导数系数放在下面循环计算
vel = 0d0
do i= 1,N+1
vel = vel + H(i)*p(i) * (i-1)
end do
end subroutine polyvel