Fortran样条插值结果异常，请求排查逻辑或公式错误

我用Fortran实现了样条插值子程序，用于超新星光变曲线的峰值时间和亮度插值，但结果存在严重偏差：

• 已知峰值亮度对应时间为53687.03535，插值得到的时间却是53639.43568
• 峰值15天后的预期亮度约18.5mag，插值结果为-5142981.63mag

观测数据：
时间（JD）：
53682.03732, 53683.04882, 53684.08633, 53687.03535, 53689.11806,
53690.06398, 53694.10385, 53695.10682, 53698.06705, 53699.09681,
53702.10265, 53706.12631, 53716.10135, 53721.06836, 53726.0874,
53730.07961, 53738.03101, 53746.03825, 53755.03675

B星等：
17.117, 17.015, 16.935, 16.838, 16.863, 16.903, 17.167, 17.25,
17.562, 17.664, 18.045, 18.583, 19.37, 19.713, 19.945, 20.141,
20.328, 20.357, 20.547

以下是实现的三个子程序：

SUBROUTINE spline(x, y, n, y1, yn, y2)
! =====================================================
! Input x and y=f(x), n (dimension of x,y), (Ordered)
! y1 and yn are the first derivatives of f in the 1st point and the n-th
! Output: array y2(n) containing second derivatives of f(x_i)
! =====================================================
IMPLICIT NONE
INTEGER:: n, i, j
INTEGER, PARAMETER:: n_max = 500
REAL*8, INTENT(in):: x(n), y(n), y1, yn
REAL*8, INTENT(out):: y2(n)
REAL*8:: p, qn, sig, un, u(n)

IF (y1 > .99e30) THEN    ! natural spline conditions
  y2(1) = 0
  u(1) = 0
ELSE
  y2(1) = -0.5
  u(1) = (3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-y1)
END IF

DO i = 2, n-1                            ! tridiag. decomposition
  sig = (x(i)-(i-1))/(x(i+1)-x(i-1))
  p = sig*y2(i-1)+2.
  y2(i) = (sig-1.)/p
  u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
END DO

IF (yn > .99e30) THEN   ! natural spline conditions
  qn = 0
  un = 0
ELSE 
  qn = 0.5
  un=(3./(x(n)-x(n-1)))*(yn-(y(n)-y(n-1))/(x(n)-x(n-1)))
END IF

y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)

DO j = n-1, 1, -1          !  backwards substitution tri-diagonale
  y2(j) = y2(j)*y2(j+1)+u(j)
END DO
RETURN

END SUBROUTINE spline

SUBROUTINE splint(x_in, y_in, spline_res, n, x_0, y_final)
! =====================================================
! Subroutine that does the actual interpolation
! Input arrays of x_in and y_in=f(x), spline_res is the result of 
! the 'spline' subroutine, x_0 is the corresponding value we are looking for
! i.e. (time_at_max in hubble), y_final is the output result
! =====================================================
IMPLICIT NONE
 INTEGER:: n, k, k_low, k_high
 REAL*8, INTENT(in):: x_in(n), y_in(n), spline_res(n), x_0
 REAL*8, INTENT(out):: y_final
 REAL*8:: a, b, h
 
 k_low = 1
 k_high = n
 
99 IF (k_high - k_low > 1) THEN
   k = (k_high + k_low) / 2
   IF (x_in(k) > x_0) THEN
     k_high = k
   ELSE
     k_low = k
   END IF
GOTO 99
ENDIF

h = x_in(k_high) - x_in(k_low) 
IF (h == 0) STOP "Bad x_in input"
a = (x_in(k_high)-x_0)/h
b = (x_0 - x_in(k_low))/h
y_final = a*y_in(k_low)+b*y_in(k_high)+((a**3-a)*spline_res(k_low)+(b**3-b)*spline_res(k_high))*(h**2)/6

RETURN 
 
END SUBROUTINE splint

SUBROUTINE spline_interp(x, y, n, x0, y_out)
! =====================================================
! Simply merging spline and splint in one subroutine
! input x and y and get y_out at x0
! =====================================================
IMPLICIT NONE
  INTEGER::n
  REAL*8, INTENT(in):: x(n), y(n), x0
  REAL*8, INTENT(out):: y_out
  REAL*8:: y1, yn, res(n)
  
! natural conditions attempt, change if not working well

y1 = 0.5
yn = 0.5
  
CALL spline(x, y, n, y1, yn, res)

CALL splint(x, y, res, n, x0, y_out)

END SUBROUTINE spline_interp

错误排查与修复 #

1. Spline子程序核心公式错误 #

在spline子程序的DO循环中，sig的计算存在致命错误：

sig = (x(i)-(i-1))/(x(i+1)-x(i-1))

这里错误地将数组索引i-1当成了前一个x值，正确的公式应该是用前一个x数组元素：

sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))

这个错误会导致sig的值完全偏离预期，进而使整个二阶导数数组y2的计算彻底失效，最终插值结果出现天文数字级的错误。

2. 边界条件设置错误 #

在spline_interp中，你设置了y1=0.5和yn=0.5，这不符合自然样条的边界条件（自然样条要求端点二阶导数为0，对应传入y1>0.99e30的标记值）。如果需要自然样条，应该修改为：

y1 = 1.0e31  ! 触发自然样条边界条件
yn = 1.0e31

如果确实需要指定端点一阶导数，需要确保输入的导数数值符合光变曲线的实际趋势（比如星等随时间的变化率），0.5的导数对于星等-时间曲线来说不合理。

3. 插值方向逻辑错误 #

你提到要寻找峰值亮度对应的时间，但当前代码是给定时间x求星等y。峰值亮度对应星等的最小值（16.838），要找对应的时间，需要做反插值：将星等作为x数组，时间作为y数组，然后插值星等=16.838对应的时间。

如果直接用现有代码传入时间求星等，得到的是该时间点的星等，而不是反向的时间值。

修复后的代码示例 #

修正后的spline子程序 #

SUBROUTINE spline(x, y, n, y1, yn, y2)
! =====================================================
! Input x and y=f(x), n (dimension of x,y), (Ordered)
! y1 and yn are the first derivatives of f in the 1st point and the n-th
! Output: array y2(n) containing second derivatives of f(x_i)
! =====================================================
IMPLICIT NONE
INTEGER:: n, i, j
REAL*8, INTENT(in):: x(n), y(n), y1, yn
REAL*8, INTENT(out):: y2(n)
REAL*8:: p, qn, sig, un, u(n)

IF (y1 > .99e30) THEN    ! natural spline conditions
  y2(1) = 0.0d0
  u(1) = 0.0d0
ELSE
  y2(1) = -0.5d0
  u(1) = (3.0d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-y1)
END IF

DO i = 2, n-1                            ! tridiag. decomposition
  sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))  ! 修正此处
  p = sig*y2(i-1) + 2.0d0
  y2(i) = (sig - 1.0d0)/p
  u(i) = (6.0d0*((y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1)) - sig*u(i-1))/p
END DO

IF (yn > .99e30) THEN   ! natural spline conditions
  qn = 0.0d0
  un = 0.0d0
ELSE 
  qn = 0.5d0
  un = (3.0d0/(x(n)-x(n-1)))*(yn - (y(n)-y(n-1))/(x(n)-x(n-1)))
END IF

y2(n) = (un - qn*u(n-1))/(qn*y2(n-1)+1.0d0)

DO j = n-1, 1, -1          !  backwards substitution tri-diagonale
  y2(j) = y2(j)*y2(j+1) + u(j)
END DO
RETURN

END SUBROUTINE spline

修正后的spline_interp子程序（自然样条） #

SUBROUTINE spline_interp(x, y, n, x0, y_out)
! =====================================================
! Simply merging spline and splint in one subroutine
! input x and y and get y_out at x0
! =====================================================
IMPLICIT NONE
  INTEGER::n
  REAL*8, INTENT(in):: x(n), y(n), x0
  REAL*8, INTENT(out):: y_out
  REAL*8:: y1, yn, res(n)
  
! 使用自然样条边界条件
y1 = 1.0e31
yn = 1.0e31
  
CALL spline(x, y, n, y1, yn, res)
CALL splint(x, y, res, n, x0, y_out)

END SUBROUTINE spline_interp

反插值求峰值时间的调用示例 #

如果要找星等=16.838对应的时间，需要交换x和y数组：

! 假设time数组存储JD时间，mag数组存储B星等
REAL*8:: time(19), mag(19), peak_time, mag_at_15days
! 初始化time和mag数组...

! 反插值：用mag作为x，time作为y，求mag=16.838对应的time
CALL spline_interp(mag, time, 19, 16.838, peak_time)

! 求峰值15天后的亮度：峰值时间是53687.03535，15天后是53687.03535+15=53702.03535
CALL spline_interp(time, mag, 19, 53702.03535, mag_at_15days)

内容的提问来源于stack exchange，提问作者Marco Leonardi